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NON-LINEAR  DYNAMIC  ANALYSES  OF  GEARED  SYSTEMS 

Rajendra  Singh,  Donald  R.  Houser,  and  Ahmet  Kahranian 

Ohio  State  University 
Department  of  Mechanical  Engineering 
Columbus,  Ohio  43210 

Under  the  driving  conditions,  a  typical  geared  system  may  be  subjected  to  large 
dynamic  loads.  Also,  the  vibration  level  of  the  geared  system  is  directly  related  to  the 
noise  radiated  from  the  gear  box.  Accordingly,  a  good  understanding  of  the  steady 
state  dynamic  behavior  of  the  system  is  required  in  order  to  design  reliable  and  quiet 
transmissions.  It  is  the  main  focus  of  this  study  with  the  scope  limited  to  a  system 
containing  a  spur  gear  pair  with  backlash  and  periodically  time-varying  mesh  stiffness, 
and  rolling  element  bearings  with  clearance  type  non-linearities.  The  internal  static 
transmission  error  excitation  at  the  gear  mesh,  which  is  of  importance  from  high 
frequency  noise  and  vibration  control  view  point,  is  considered  in  the  fonnulation  in 
sinusoidal  or  periodic  form. 

A  dynamic  finite  element  model  of  the  linear  time-invariant  (LT1)  system  is 
developed.  Effects  of  several  system  parameters,  such  as  torsional  and  transverse 
flexibilities  of  the  shafts  and  prime  mover/load  inertias,  on  free  and  forced  vibration 
characteristics  are  investigated.  Several  reduced  order  LTI  models  are  developed  and 
validated  by  comparing  their  eigen  solutions  with  the  finite  element  model  results. 
Using  the  reduced  order  formulations,  a  three-degree  of  freedom  dynamic  model  is 
developed  which  includes  non-linearities  associated  with  radial  clearances  in  the  radial 
rolling  element  bearings,  backlash  between  a  spur  gear  pair  and  periodically  varying 
gear  meshing  stiffness.  As  a  limiting  case,  a  single  degree  of  freedom  model  of  the 
spur  gear  pair  with  backlash  is  considered  and  mathematical  conditions  for  tootli 
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separation  and  back  collision  are  defined.  Both  digital  simulation  technique  and 
analytical  methods  such  as  method  of  harmonic  balance  and  the  method  of  multiple 
scales  have  been  used  to  develop  the  steady  state  frequency  response  characteristics  for 
various  non-linear  and/or  time-varying  cases.  Difficulties  associated  with  the 
detennination  of  the  multiple  solutions  at  a  given  frequency  in  the  digital  simulation 
technique  have  been  resolved.  The  proposed  formulation  has  been  validated  by 
comparing  the  predictions  with  the  results  of  two  benchmark  experiments  reported  in 
the  literature.  Several  key  system  parameters  such  as  mean  load  and  damping  ratio  are 
identified  and  their  effects  on  the  non-linear  frequency  response  are  evaluated 
quantitatively.  Other  fundamental  issues  such  as  the  dynamic  coupling  between  non¬ 
linear  modes,  dynamic  interactions  between  component  non-linearities  and  time-varying 
mesh  stiffness,  and  the  existence  of  subhannonic  and  chaotic  solutions  including  routes 
to  chaos  have  also  been  examined  in  depth. 
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CHAPTER  I 


INTRODUCTION 

1.1.  PROBLEM  FORMULATION 

Dynamic  analysis  of  geared  systems  is  an  essential  step  in  design  due  to  two 
reasons.  First,  under  the  driving  conditions,  a  typical  geared  system  is  subject  to 
dynamic  forces  which  can  be  large.  Therefore,  the  prediction  of  dynamic  loads, 
motions  or  stresses  is  needed  in  developing  reliable  gear  trains.  Second,  the  vibration 
level  of  the  geared  system  is  directly  related  to  the  noise  radiated  from  the  gear  box  An 
attempt  in  designing  quiet  gears  requires  a  good  understanding  of  the  dynamic  behavior 
of  the  system  and  the  gear  mesh  source.  Accordingly,  the  main  objective  of  this  study 
is  to  develop  accurate  mathematical  models  of  a  generic  geared  rotor-bearing  system 
shown  in  Figure  1.1a.  Of  interest  here  is  to  investigate  several  key  modelling  issues 
which  have  not  been  addressed  in  the  literature,  such  as  system  non-linearities  and 
time-varying  mesh  stiffness. 

The  generic  geared  system  shown  in  Figure  1.1a  consists  of  a  single  spur  gear 
mesh  of  ratio  vg=dg2/dgl,  rolling,  element  bearings,  a  prime  mover  driving  the  system 

at  Osi  speed  and  a  typical  inertial  load.  The  system  also  includes  other  elements  such 

as  couplings  and  flywheel.  A  discrete  model  of  the  system  is  shown  in  Figure  1.1b. 
Here,  shafts  are  represented  by  discrete  translational  springs  ksJ  and  ks2,  translational 
dampers  csl  and  cs2,  and  torsional  springs  KsJ  and  K$2  The  gear  mesh  is  represented 
by  a  time-varying  mesh  stiffness  kh(t)  and  a  non-linear  displacement  function  fh 


Ig2.nig2.dg2 


(b) 


Figure  1.1.  a)  A  generic  geared  rotating  system,  b)  discrete  model  of  geared  rotating 


which  includes  gear  backlash.  Further,  linear  time-invariant  (LTI)  mesh  damping  cb  is 

considered  here.  The  rolling  element  bearings  are  defined  by  a  time-invariant  radial 
stiffness  kb  subject  to  a  non-linear  displacement  function  fb,  and  an  LTI  damping 
coefficient  cb.  The  prime  mover  and  load  are  modeled  as  purely  torsional  elements  of 
inertias  Ip  and  IL,  respectively.  The  mean  rotational  speeds  Q8j  and  ftS2  and  the 

geometric  end  conditions  are  such  that  gyroscopic  effects  are  not  seen. 

The  generalized  displacement  vector  {q(t)},  associated  with  the  inertia  elements, 

consists  of  angular  displacements  0  and  transverse  displacements  x  and  y.  The 
governing  equation  of  motion  for  the  non-linear,  time-  varying  multi-degree  of  freedom 
model  can  be  given  in  the  general  form  as 

[M]{q"(  t )}  +  [c]{q'(t )}  +  [K(i)]{f  (q(t ))}  =  {F(t )}  (1.1) 

where  [m]  is  the  time-invariant  mass  matrix  and  {q(t)}  is  the  displacement  vector. 
Here,  damping  matrix  [cj  is  assumed  to  be  LTI  type,  as  the  effect  of  the  tooth 
separation  and  time-varying  mesh  properties  on  mesh  damping  are  considered 

negligible;  validity  of  this  assumption  will  be  examined  later.  The  stiffness  matrix 
[K(0]  is  considered  to  be  time-varying,  given  by  a  periodically  time  varying  matrix 

[K(t)]  =  [K(t  +  2n/f2j,)]  where  Qh  is  the  fundamental  gear  mesh  frequency.  The 

non-linear  displacement  vector  {f(q(t))}  includes  the  radial  clearances  in  bearings  and 
the  gear  backlash,  and  the  forcing  vector  |F(t)}  consists  of  both  external  excitations 

due  to  torque  fluctuations,  mass  unbalances  and  geometric  eccentricities,  and  an 
internal  static  transmission  error  excitation. 
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1.2.  OBJECTIVES 

Specific  objectives  of  this  study  are  given  as  follows;  each  chapter,  written  in  the 
journal  paper  style,  deals  with  one  major  objective. 

First,  a  dynamic  finite  element  model  of  the  linear  time-invariant  (LTI)  system 
given  in  Figure  1.1a  is  developed.  Effects  of  several  system  parameters,  such  as 
torsional  and  transverse  flexibilities  of  the  shafts  and  prime  mover/load  inertias,  on  free 
and  forced  vibration  characteristics  are  investigated.  Three  different  reduced  order  LTI 
models  will  be  derived  and  the  conditions  under  which  these  are  valid  will  be 
determined  by  comparing  the  eigen  solutions  with  the  finite  element  model  results. 
Development  and  verification  of  such  a  reduced  order  (with  a  very  few  degrees  of 
freedom)  linear  model  is  an  essential  step  before  the  non  linear  dynamic  behavior  is 
analyzed  [Chapter  I], 

Second,  non-linear  frequency  response  characteristics  of  a  spur  gear  pair  with 
backlash  are  examined  for  both  external  and  internal  excitations.  The  internal  excitation 
is  of  importance  from  the  high  frequency  noise  and  vibration  control  viewpoint  and  it 
represents  the  overall  kinematic  or  static  transmission  error.  Such  problems  may  be 
significantly  different  from  the  rattle  problems  associated  with  external,  low  frequency 
torque  excitation.  Two  solution  methods,  namely  the  digital  simulation  technique  and 
the  method  of  harmonic  balance  have  been  used  to  develop  the  steady  state  solutions 
for  the  internal  sinusoidal  excitation.  Difficulties  associated  with  the  determination  of 
the  multiple  solutions  at  a  given  frequency  in  the  digital  simulation  technique  have  been 
resolved  as  one  must  search  the  entire  initial  conditions  map.  Such  solution?  and  the 
transition  frequencies  for  various  impact  situations  are  found  by  »he  method  of 
harmonic  balance.  Further,  the  principle  of  superposition  can  be  employed  to  analyze 
the  periodic  transmission  error  excitation  and/or  combined  excitation  problems 
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provided  the  excitation  frequencies  are  sufficiently  apart  from  each  other.  Predictions 
are  compared  with  the  limited  experimental  data  available  in  the  literature  [Chapter  II]. 

Third,  non-linear  frequency  response  characteristics  of  a  geared  rotor-bearing 
system  are  examined.  A  three  degree  of  freedom  dynamic  model  is  developed  which 
includes  non-linearities  associated  with  radial  clearances  in  the  radial  rolling  element 
bearings  and  backlash  between  a  spur  gear  pair;  time-invariant  gear  meshing  stiffness  is 
assumed.  The  bearing  non-linear  stiffness  function  is  approximated  for  convenience 
sake  by  a  simple  model  which  is  identical  to  that  used  for  the  gear  mesh.  This 
approximate  bearing  model  has  been  verified  by  comparing  steady  state  frequency 
spectra.  The  applicability  of  both  analytical  and  numerical  solution  techniques  to  the 
multi  degree  of  freedom  non-linear  problem  is  investigated.  Proposed  theory  is 
validated  by  comparing  the  results  with  available  experimental  data.  Several  key 
issues  such  as  non-lincar  modal  interactions  and  differences  between  internal  static 
transmission  error  excitation  and  external  torque  excitation  are  discussed.  Additionally, 
parametric  studies  are  performed  to  understand  the  effect  of  system  parameters  such  as 
bearing  stiffness  to  gear  mesh  stiffness  ratio,  alternating  to  mean  force  ratio  and  radial 
bearing  preload  to  mean  force  ratio  on  the  non-linear  dynamic  behavior.  A  criterion 
used  to  classify  the  steady  state  solutions  is  presented  and  the  conditions  for  chaotic, 
quasi-periodic  and  subharmonic  steady  state  solutions  are  determined.  Two  typical 
routes  to  chaos  observed  in  this  geared  system  are  also  identified  [Chapter  ID]. 

Fourth,  this  study  extends  the  non-linear  single  degree  of  freedom  spur  gear  pair 
model  of  Chapter  II  and  multi-degree  of  free  dom  geared  rotor-bearing  system  model  of 
Chapter III  by  including  lime-varying  mesh  stiffness  k^tt),  and  investigates  the  effect 

of  k h (t )  on  the  frequency  response  of  lightly  and  heavily  loaded  geared  systems. 
Interactions  between  mesh  stiffness  variation  and  system  non-linearities  associated  with 
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gear  backlash  and  radial  clearances  in  rolling  element  bearings  are  also  considered. 
Resonances  of  the  corresponding  linear  time-varying  (LTV)  system  associated  with  the 
parametric  and  external  excitations  are  identified  using  the  method  of  multiple  scales. 
Theoretical  results  are  validated  by  available  experimental  results  (Chapter  IV). 

1.3.  DEVELOPMENT  OF  LINEAR  TIME-INVARIANT  MODELS 
1.3.1.  Literature  Review 

The  study  of  geared  rotor  dynamics  requires  that  the  coupling  between  torsional 
and  transverse  vibrations  be  included  in  the  model.  Although  several  modeling  and 
solution  techniques  such  as  lumped  mass  models  and  the  use  of  the  transfer  matrix 
method  have  been  applied  to  rotor  dynamic  problems,  the  Finite  element  method  (FEM) 
seems  to  be  a  highly  efficient  and  accurate  method  for  linear  modeling.  In  one  of  the 
early  examples  of  FEM  applied  to  a  single  rotor,  Nelson  and  McVaugh  [1]  used  a 
Rayleigh  beam  finite  element  including  the  effects  of  translational  and  rotary  inertia, 
gyroscopic  moments  and  axial  load.  Zorzi  and  Nelson  [2]  extended  this  by  including 
internal  damping  Later,  Nelson  [3]  developed  a  Timoshenko  beam  by  adding  shear 
deformation  to  the  Rayleigh  beam  theory.  This  model  was  further  extended  by 
Ozguven  and  Ozkan  (4)  to  include  all  possible  effects  such  as  transverse  and  rotary 
inertia,  gyroscopic  moments,  axial  load,  internal  hysteretic  and  viscous  damping  and 
shear  deformations  in  a  single  model.  However,  none  of  the  rotor  dynamics  models 
described  above  can  handle  geared  rotor  systems,  although  they  are  capable  of 
determining  the  dynamic  behavior  of  rotors  which  consist  of  shafts  supported  at  several 
points  and  carrying  rigid  disks  at  several  locations. 

Gear  dynamics  studies,  on  the  other  hand,  have  usually  neglected  the  lateral 
vibrations  of  the  shafts  and  bearings,  and  have  typically  represented  the  system  with  a 


linear  torsional  model.  Although  neglecting  transverse  vibrations  might  be  a  good 
approximation  for  systems  having  stiff  shafts,  it  has  been  observed  experimentally  [5] 
that  the  dynamic  coupling  between  the  transverse  and  torsional  vibrations  due  to  the 
gear  mesh  affects  the  system  behavior  considerably  when  the  shafts  are  compliant. 
This  fact  directed  the  attention  of  investigators  to  the  inclusion  of  transverse  vibrations 
of  the  shafts  and  the  bearings  in  mathematical  models.  Lund  [6]  developed  influence 
coefficients  at  each  gear  mesh  for  both  torsional  and  lateral  vibrations,  obtained  critical 
speeds  and  a  forced  vibration  response.  Hamad  and  Seireg  [7]  studied  the  whiiling  cf 
geared  rotor  systems  supported  on  hydrodynamic  bearings;  torsional  vibrations  were 
not  considered  in  this  model  and  the  shaft  was  assumed  to  be  rigid.  Iida,  et  al.  [8] 
considered  the  same  problem  by  taking  one  of  the  shafts  to  be  rigid  and  neglecting  the 
compliance  of  the  gear  mesh,  and  obtained  a  three  degree  of  freedom  model  to 
determine  the  first  three  vibration  modes  and  the  forced  vibration  response  due  to  the 
unbalance  and  the  geometric  eccentricity  of  one  of  the  gears.  Later,  Iida,  et  al. 
[9,10,1 1]  applied  their  model  to  a  larger  system  which  consists  of  three  shafts  coupled 
by  two  gear  meshes.  Hagiwara,  Iida,  and  Kikuchi  [12]  developed  a  simple  model  that 
included  the  transverse  flexibilities  of  the  shafts  by  using  discrete  stiffness  values,  and 
studied  the  forced  response  of  geared  shafts  due  to  unbalances  and  runout  errors.  They 
included  the  damping  and  compliances  of  the  journal  bearings  and  assumed  a  constant 
mesh  stiffness.  Although  most  of  these  gear  dynamics  studies  have  discussed  several 
aspects  of  the  problem,  none  of  them  has  been  able  to  represent  a  geared  rotor-bearing 
system  completely  since  almost  all  of  them  have  used  one  or  more  simplifying 
assumptions  such  as  rigid  shafts,  rigid  bearings,  rigid  gear  mesh  etc.  which  may  not  be 
applicable  to  a  real  system.  In  addition,  these  studies  proposed  lower  order  lumped 
mass  linear  models. 
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Neriya,  et  al.  [13]  employed  a  dynainic  finite  element  model  which  eliminates 
many  of  the  simplifying  assumptions.  They  found  the  forced  vibration  response 
of  the  system  at  the  shaft  frequency,  excited  by  mass  unbalances  and  runout  errors  of 
the  gears  by  using  the  modal  summation.  But  they  did  not  consider  the  high 
frequency,  internal,  static  transmission  error  excitation  which  has  the  major  role  in 
noise  generation.  An  extensive  survey  of  linear  mathematical  models  used  in  gear 
dynamics  analyses  is  given  in  a  recent  paper  by  Ozguven  and  Houser  [14], 

1.3.2.  Mathematical  Model 

Here,  we  assume  linear  bearings  and  no  tooth  separations  with  time-invariant 
mesh  stiffness,  i.e.  [k]  *  [lC(t )]  and  {f(q(t ))}  =  {q(t)}.  Then  the  corresponding  LTI 

form  of  equation  ( 1 . 1 )  is: 

[M]{q"(i )}  +  [c]{q'(t )}  +  [K]{q(t)}  =  {F(i>)  (12) 

Since  many  investigators  have  modeled  the  generic  system  shown  in  Figure  1.1a  as  a 
single  degree  of  freedom  model,  our  analysis  uses  it  as  a  reference  model  to  transform 
the  governing  equation  Into  the  dimensionless  form.  Focusing  only  on  the  gear  pair  as 
shown  in  Figure  1.2,  the  equation  of  motion  of  the  semi-definite  system  is  given  in 
terms  of  the  relative  translational  displacement  u(t)  =  (dgjGgi  -  dg2©g2)  /  2 

d2u(t)  du(t)  ,  x 

mci— rri  +  ch-^i  +  khu(0  =  F(t) 
dt  dt 

where  is  the  equivalent  gear  pair  mass  defined  as 
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1 


(1.3b,c) 


mcl  = 


j2  .2 

dgl  ,  dg2 


4Ig,  41 


g2 


First,  we  establish  the  dimensionless  time  t  as  t  =  o>nt .  Second,  we  use  the  base  circle 
diameter  of  the  pinion  dgj  and  the  equivalent  mass  mcj  as  characteristic  length  and 
mass  parameters,  respectively,  to  obtain  the  governing  equations  in  the  dimensionless 
form  as 


[M]{q(t)J  +  [C){q(t)}  +  [K]{q(t)}  =  {Fm}  +  {Fah(t)}  +  {F,T(l)} .  (1.4) 

where  an  overdot  means  derivative  with  respect  to  time  t,  and  the  dimensionless  forcing 
vector  consists  of  a  mean  force  vector  |Fm )  and  two  time-varying  components:  a)  high 
frequency  excitation  due  to  the  kinematic  gear  transmission  error  IF^ft)}  and.  b)  other 
excitations  due  to  mass  unbalance  Ugj,  geometric  eccentricities  (run-out  errors)  £gj  and 
p  ime  mover  and  load  torque  pulsations  Tgia(t),  typically  at  low  frequencies,  are 
combined  into  a  single  term  (FaT(t)) . 

1.3.3.  Gear  Mesh  Formulation 

The  gear  mesh  is  represented  by  a  pair  of  rigid  disks  connected  by  a  translational, 
viscously  damped  spring  along  the  pressure  line  which  is  tangent  to  the  base  circles  of 
the  gears  as  shown  in  Figure  1.3.  By  choosing  the  Y  axis  on  the  pressure  line  and  the 
X  axis  perpendicular  to  the  pressure  line,  the  transverse  vibrations  in  the  X  direction 
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are  uncoupled  from  both  the  torsional  vibrations  and  transverse  vibrations  in  the  Y 
direction.  The  dynamic  mesh  force  Fby  at  the  mesh  point  in  the  Y  direction  can  be 

written  in  terms  of  the  symbols  given  in  Figure  1 .3  from  the  rigid  body  dynamics 
Fhvlt) = ch(ygi + ^0ti + egiftsi  cos0ti  -  yg2  -  -^2 

dei 

-eg2^s2cos0t2  -eNglQslcos(Ngl0ti))+kh(ygi  +  — 9,i 
+£gisin0tl  -  yg2 - jp-0,2  “£g2sin0t2  -  ej  sin(Ngl0,i ))  U-5) 

Here,  only  the  fundamental  component  of  the  static  transmission  error  e(t)  is 
considered,  i.e.  e(t)=  esin(NglQsit)  where  Ng)  is  the  number  of  teeth  of  gear  gl. 
Dynamic  mesh  force  Fhy  also  introduces  moments  ’f'hgl  and  Thg2  about 
instantaneous  centers  of  the  gears: 

Thgi(t)  =  +  £gi  cos9tl]  0  -6a) 

Thg2(i)  =  +  eg2  COS0,2  j  ( 1  -6b) 

The  total  angular  rotation  0t,  of  the  i-th  gear  is  ©tj (t )  =  flsjt  +  0gi(t)  ,  i=l,2.  The 
displacement  vector  can  be  decomposed  into  two  parts  as  {q}  =  {%}  +  {^h}  where 

{qs}  =  [MO?l.0l . . ^n.yD.0n].  j=lton,  j*gl,j*g2  (excluding  the 

degrees  of  freedom  of  the  nodes  where  the  gear  and  pinion  are  mounted)  and 
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{qh}  =  [xgi,ygi,eg|,Xg2,yg2,0g2],  i=l ,2.  The  coupling  effect  due  to  gear  mesh  is 
only  seen  in  the  terms  governed  by  {qh}  and  the  mesh  stiffness  [Kh],  which 
represents  this  dynamic  coupling  due  to  gear  mesh  is  obtained  from  equations  (1.5)  and 
(1.6)  as 


[Kb]=  k. 


0 

0 

0 

0 

0 

0 


0 

1 

£gi 

2 

0 

-1 


Jg2 


0 

dgl 

2 

4 

0 

dp  1 
2 

dgidg2 


0 

0 

0 

0 

0 

0 


0 

-1 

.£gi 

2 

0 

1 

dg2 


0 

dg2 

2 

dgldg2 


0 

Jg2 


*g2 


(17a) 


Since  mesh  viscous  damping  is  proportional  to  the  mesh  stiffness, 


(1.7b) 


In  the  overall  force  vector,  non-zero  terms  correspond  to  displacements  {qh}  and  are 
given  in  the  following  form: 
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{Fah(‘)}  +  {FaT(‘)}  =  | 


(1.8a) 


UgA,  sinfis|t  4-  F(t)j 

uglQsicosnsit 

-£Zsl2lE  cosiigit  +  F,(i)  +  Te„(t) 

dgl  2 

Ug2«?2sinns2t-^(t) 

ug2^?2  cosQs2t 

"1~£g2ns2  cosfts2t - Fj  (t)  -  Tg2a(t) 

a0i  1 


Fi(t)  =  Cj1(£g2t2s2  cosOs2t  £g|^si  cosaslt  4-  e  Ngjlis|  cos(Ngjf2sjt )) 

+kb(eg2  sinfts2t  -  Eg]  sinQslt  +  esin(Nglftslt ))  ( 1  8b) 

where  Tg]m  is  the  mean  input  torque. 

1.3.4.  Finite  Element  Formulation  and  Eigen-Value  Problem 

The  finite  element  method  has  been  employed  in  obtaining  stiffness  and  mass 
matrices  of  the  system  of  equations  ( 1 .4).  The  shafts  are  discretized  and  five  degrees  of 
freedom  are  defined  at  each  node,  only  the  axial  motion  being  excluded.  The  stiffness 
and  mass  matrices  of  each  finite  rotor  element  are  derived  by  using  the  variational 
principle  [3,4,15].  The  system  overall  matrices  are  obtained  by  combining  element 
matrices,  dynamic  coupling  matrices  due  to  gear  mesh  defined  by  equation  (1.7), 
assumed  bearing  damping  values  and  stiffnesses,  and  the  inertias  of  the  gears  and  other 
lumped  inertia  elements.  Here,  [M]  is  diagonal  and  positive  definite  and  [K]  is 
symmetric  and  positive  semi-definite.  The  equation  describing  the  undamped  free 
vibration  of  the  system  is  obtained  from  equation  (1 .4)  as 
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[M]{q(t)}  +  [K]{q(t)}  =  {0}. 


(19) 


The  assumed  form  of  the  solution  (q(t)}  =  {qa}cos(o)t  +  $)  is  substituted  into  equation 
(1.9)  to  obtain  the  standard  eigen-value  problem: 

[K]{qa}  =  CU2[M]{qa}.  (1.10) 

The  solution  of  equation  (1.10)  yield  the  eigen-vectors  or  modes  (\j/r)  and  associated 
eigen-values  or  natural  frequencies  tor.  Here,  a  torsional  rigid  body  mode  at  co=0 
exists  since  [K]  is  semi-definite. 

1.3.5.  Forced  Vibration  Response 

The  excitation  given  by  equation  (1.8)  (for  no  torque  pulsations,  i.e. 
Tgja(t)  =  Tg2a(t)  =  0)  is  the  sum  of  three  sinusoidal  terms  at  frequencies  Qs|,  and 

gear  mesh  frequency  Ngj£2si 

(p(0}=  £  {Fg,}^«s,<+4’siHM^NsAi'+'t’b)  on) 

i=l 


The  steady  state  displacement  response  of  the  system  due  to  this  excitation  is  assumed 
to  be 


{q(*)}  =  I  [Pp  ]{Fg,}r,in(fts,t  +  4>S1 )  +  [Ph]{Fh}sm(Nglnslt  +  <J>h)  (1.12) 

i  =  l 
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where  (figj]  i=l,  2  and  [p^]  are  the  dynamic  compliance  matrices  in  the  frequency 
domain  corresponding  to  the  excitation  frequencies,  ns|,fts2and  NgjQsl, 
respectively. 


(¥r)l¥r)T 
(ror  -  ^si) + 


i  =  1,  2, 


(1.13a) 


l¥r)(Vrl 


r=l 


(0r  ^Ng|f2sjj  j 2j^rNgjQSjtur 


(1.13b) 


{ q/r )  represents  the  r-th  mass  matrix -normalized  modal  vector,  n  is  the  total  number  of 
degrees  of  freedom  of  the  system,  j  =  V-l ,  and  £r  is  the  damping  ratio  for  r-th  mode. 
See  Appendix  A  for  the  computer  code  GRD  which  uses  the  theory  given  here. 

1.4.  PARAMETRIC  STUDIES 
1.4.1.  Modes  of  Interest 

The  generic  system  shown  in  Figure  1.1  is  modeled  by  FEM  to  examine  the 
natural  modes  of  a  general  linear  geared  rotor-bearing  system  and  to  study  the  effects  of 
several  system  parameters  on  the  dynamics  of  the  system.  Two  numerical  data  sets  as 
listed  in  Table  1.1  are  used.  Predicted  natural  frequencies  and  modes  for  Set  A  are 
presented  in  Table  1.2.  The  response  of  the  system  to  e(t)  is  also  computed;  Figures 

1.4,  1.5  and  1.6  display  the  response  in  the  Y  and  torsional  directions  at  the  pinion 
location  and  the  dynamic  load  to  static  load  ratio  at  the  mesh  point  dgiF^y  /  2Tmgl , 

respectively.  The  system  has  no  peak  responses  at  the  modes  corresponding  to  motion 
in  the  X  direction  since  the  excitation  is  applied  in  the  Y  (pressure  line)  direction,  and 


1(> 


Table  1.1  Numerical  data  sets  of  the  system  used  for  calculations. 


Parameters 

Set  At 

Set  B 

Igl,Ig2  (kg-m2) 

0.0018,  0.0018 

0.009/,  0.0097 

mgl,mg2  (kg) 

1.84,  1.84 

3.45,  3.45 

dg[,  dg2  (m) 

0.089,  0.089 

0.135,  0.135 

Ngl 

28 

30 

kh 

1.0x10® 

1.0x10® 

kb  (N/m) 

variable 

rigid 

LslL-LslR  <m) 

0.127,  0.127 

variable 

Ls2L>Ls2R  (m) 

0.127,  0.127 

variable 

KltK2  (N-m/rad) 

— 

variable 

Ip.  IL  (kg-m2) 

— 

variable 

c4l,0’^s2,o 

0.037,  0.037 

0.04,  0.04 

^sl  ,i'  ^s2,i 

0.,  0.01 

— 

«  (m) 

9.3xl0‘6 

-- 

t  NASA  Lewis  Research  Center  gear  test  rig. 
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Tabic  1.2.  First  10  natural  frequencies  for  set  A  of  Table  1.1  (kfc/kb=10). 


Modal  Modes  of  Natural  Frequency  to,  Natural  Mode 

Index  Interest  Hz.  Description 


0 

0 

torsional  rigid  body 

1 

Vi 

581 

first  transverse-torsional  coupled 

2 

687 

X  direction,  transverse  (driving  shaft) 

3 

Vii 

689 

Y  direction,  transverse 

4 

691 

X  direction,  transverse  (driven  shaft) 

5 

Vuj 

2524 

second  transverse-torsional  coupled 

6 

3387 

Y  direction,  transverse 

7 

3387 

X  direction,  transverse  (driving  shaft) 

8 

3421 

X  direction,  transverse  (driven  shaft) 

9 

3421 

Y  direction,  transverse 

18 


in  Y -direction  (m) 


Figure  1.4.  Forced  response  of  the  system  with  dataset  A  (Table  1.1)  to  the 
displacement  excitation  in  the  direction  of  pressure  line  (at  pinion 
location)  for  four  different  bearing  stiffnesses  kb. 
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torsional  response  (Rad.) 


|H>- 


dynamic  to  static  load  ratio 


2.0 

1.8 

1.6 

1.4 

1.2 

1.0 

excitation  frequency  (Hz) 

Figure  1.6  Dynamic  to  static  load  ratio  for  dataset  A  (Table  1.1)  due  to  the  static 
transmission  error  excitation  for  four  different  beaiing  stiffnesses  k^. 


500  1000  1500  2000  2500  3000 


the  vibration  in  the  X  direction  is  dynamically  uncoupled  from  the  vibrations  in  the  Y 
and  torsional  motions,  as  described  in  the  previous  section.  Therefore,  the  natural 
modes  corresponding  to  motions  in  the  X  direction  can  be  eliminated  when  only  e(t) 

excites  the  system.  Accordingly,  the  following  three  typical  modes  are  of  special 
importance: 

I.  First  Transverse-Torsional  Coupled  Mode  wp.  This  is  the  second  mode  listed 
in  Table  1 .2  at  coj  (581  Hz  for  the  system  considered).  This  mode  corresponds  to  the 
first  peak  in  Figures  1.4  and  1.5  and  its  schematic  shape  is  shown  in  Figure  1.7a. 
Here,  shafts  move  in  opposite  directions,  gears  vibrate  in  opposite  directions  also;  but 
transverse  and  torsional  vibrations  combine  to  yield  small  relative  motion  at  the  gear 
mesh  point.  Therefore,  dynamic  loads  at  the  mesh  point  are  not  large,  resulting  in  no 
peak  in  Figure  1.6  at  tOj  while  Figures  1.4  and  1.5  have  peaks  governed  by  this  mode. 

II.  Purely  Transverse  Mode  At  this  mode  with  natural  frequency  (0]j  (689 
Hz  in  Table  1.2),  there  is  no  torsional  vibration  and  both  shafts  vibrate  in  phase  in  the 
pressure  line  direction  as  shown  in  Figure  1.7b  schematically.  The  relative 
displacement  at  mesh  point  is  zero  since  the  gear  ratio  vg=l.  Therefore,  e(t)  cannot 

excite  this  mode  as  no  peaks  are  observed  in  Figures  1.4,  1.5  and  1.6  at  frequency  cou- 

III.  Second  Transverse-Torsional  Coupled  Mode  ,YhU  The  second  and  the 
highest  peak  seen  in  Figures  1.4,  1.5  and  1.6  corresponds  to  this  mode  at  o>m  .  As 
seen  from  the  mode  shape  illustrated  in  Figure  1.7c,  both  shafts  and  gears  vibrate  in 
opposite  directions,  and  transverse  and  torsional  vibrations  are  additive  at  the  mesh 
point.  Thus  a  large  relative  displacement  at  the  mesh  point  is  obtained,  which  results  in 
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Figure  1 .7.  Typical  natural  modes  of  interest;  see  Table  1 .2  for  further  details.. 
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a  large  peak  in  Figures  1.4,  1.5  and  1.6.  This  is  the  mode  at  which  the  coupling 
between  transverse  and  torsional  vibrations  is  very  strong. 

These  three  modes  are  observed  in  all  geared  rotor  systems  and  they  play  an 
important  role  in  governing  dynamic  response  of  the  system  excited  by  e(t). 

Therefore,  any  accurate  mathematical  model  of  the  geared  rotor-bearing  systems  must 
be  able  to  predict  these  modes. 

1.4.2.  Effect  of  Bearing  Compliances 

In  most  cases,  radial  stiffness  of  a  typical  rolling  element  bearings  kfe  is  roughly 
of  the  same  order  of  magnitude  as  the  gear  mesh  stiffness  kh;  in  general  it  lies  in  the 
range  0.1kh<kb<100kh.  Therefore,  bearing  flexibility  should  be  included  in  the 

analysis.  A  parametric  study  for  data  set  A  of  Table  1.1  has  been  conducted  to 
demonstrate  the  effect  of  kb  on  the  natural  frequencies  and  the  frequency  response  of 
the  system  excited  by  e(t).  Figures  1.4,  1.5  and  1.6  show  the  response  in  the  Y  and 
torsional  directions  and  the  dynamic  mesh  load  respectively  for  k^  values  ranging  from 
0.1kh  to  10kh,  shaft  stiffnesses  are  being  kept  the  same.  An  increase  in  bearing 

stiffness  results  in  an  increase  in  both  natural  frequencies  and  the  peak  amplitudes. 
Above  the  range  of  kb,  the  bearing  becomes  very  rigid  when  compared  to  the  shaft 

compliances  and  its  effect  on  the  system  can  be  ignored. 

1.4.3.  Effect  of  Shaft  Compliance 

Data  set  B  of  Table  1.1  has  been  used  to  study  the  effect  of  shaft  compliance  on 
the  natural  modes.  The  shaft  length  is  varied  and  the  corresponding  cor  values  are 
predicted  by  FEM;  Figure  1.8  displays  this  result.  The  first  natural  frequency  (Oj  does 
not  change  considerably  with  varying  shaft  length,  whereas  (Ojj  and  (Din  are  strongly 


24 


dimensionless  natural  frequency 


dependent  on  the  shaft  length,  especially  at  smaller  lengths.  Another  observation  from 
Figure  1 .8  is  that  o>ji  and  coyj  become  very  large  for  smaller  shaft  lengths  and  clearly 
move  beyond  the  range  of  the  operational  speed  of  most  geared  rotor  systems.  In  this 
case,  the  first  coupled  transverse-torsional  mode  \|fj  can  be  assumed  to  be  uncoupled 
from  these  two  modes. 

1.4.4.  Effect  of  Load  and  Prime  Mover  Rotary  Inertias 

As  shown  in  Figure  1 . 1 ,  the  shafts  are  connected  to  a  prime  mover  and  a  load  at 
either  end  through  flexible  torsional  couplings.  The  following  parameters  need  to  be 
considered:  motor  and  load  rotary  inertias,  torsional  compliances  of  the  flexible 
couplings,  and  stiffnesses  of  the  driving  and  driven  shafts. 

First,  the  motor  and  load  are  assumed  to  be  connected  to  the  shafts  without 
considering  any  torsional  couplings  in  between.  Figure  1 .9  displays  the  variation  in 
coj,  cojj  and  coju  with  a  variation  in  the  prime  mover  inertia.  Here,  data  set  B  of  Table 
1.1  with  Lsi=0.04  m  is  used.  (Oy  and  toy]  are  not  affected  and  therefore,  the  inertias  of 

motor  and  load  can  be  disregarded  if  the  major  concern  is  to  predict  these  two  modes. 
However,  (Oj  is  strongly  dependent  on  prime  mover  and  load  inertias. 

Second,  the  load  and  prime  mover  inertias  are  fixed  (Ip=5Igl)  and  the  torsional 
springs  Kj  and  K2,  which  represent  flexible  couplings  and  shafts  are  varied.  Figure 
1.10  shows  the  variation  in  (0,.  with  changing  Kj  and  K2.  Here  (Oy  and  (i>yj  are  again 
not  affected,  as  expected.  On  the  other  hand,  o)j  is  almost  constant  (which  is  nearly 
equal  to  the  value  yielded  by  zero  prime  mover  and  load  inertias)  up  to  a  point  and,  it 
starts  increasing  with  increasing  K.  This  indicates  that  the  motor  and  load  are  isolated 
from  the  geared  rotor  system  when  the  in-between  torsional  elements  are  compliant 
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dimensionless  natural  frequency 


dimensionless  natural  frequency 


°'510’  10‘  10’  10* 
torsional  stiffness  (N.m/rad) 


Figure  1.10.  Effect  of  the  torsional  stiffness  K  of  the  transmission  elements  in 
between  the  gear  box  and  the  prime  mover  &  load  inertias  on  the  natural 
frequencies. 
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enough,  which  is  the  case  in  most  practical  systems.  Under  these  circumstances, 
motor  and  load  inertias  can  be  neglected  in  the  analysis. 

1.5.  REDUCED  ORDER  LINEAR  TIME-INVARIANT  MODELS 

In  this  section,  three  different  reduced  order  analytical  models  of  the  geared  rotor¬ 
bearing  systems  shown  in  Figure  1 . 1  will  be  developed  and  the  conditions  and  system 
parameters  at  which  these  simple  models  can  predict  the  dynamics  of  the  system 
accurately,  will  be  discussed.  The  finite  element  model  of  Section  1.3  will  be 
employed  as  a  reference  model  to  check  the  validity. 

1.5.1.  Single  Degree  of  Freedom  Torsional  Model  of  dear  Pair 

As  the  simplest  model,  a  single  degree  of  freedom  (SDOF)  model  of  geared  rotor 
systems  shown  in  Figure  1.2  is  considered.  The  shaft  and  bearing  flexibilities  and  the 
motor  and  load  inertias  are  not  considered  in  this  model.  This  model  can  only  predict  a 
single  mode  at  a  frequency  ton  as  defined  in  equation  (1.3c)  which  corresponds  to  the 
first  transverse-torsional  mode  at  toj  .  Here,  o>j-»con  when  shaft  lengths  Lsi— >0  and 
bearings  are  very  stiff.  And,  (Ou  and  o>ni  are  sufficiently  beyond  the  operational  speed 
range  and  the  variation  in  (Oj  is  assumed  to  be  negligible  with  shaft  length  as  shown  in 
Figure  1 .8.  Therefore,  in  some  cases,  the  SDOF  model  of  Figure  1.2  can  be  used  to 
represent  the  system  provided  these  conditions  are  met.  For  instance,  for  data  set  B  of 
Table  1 . 1  with  1^=5  cm,  a  SDOF  model  can  be  utilized  up  to  an  operational  rotational 
speed  of  6000  rpm  which  corresponds  to  the  excitation  frequency  at  3000  Hz  for 
Ngt=30  teeth.  As  it  is  seen  in  Figure  1.8,  only  the  first  mode  is  observed  and  its 
variation  is  not  significant  within  0£f2<3000  Hz  and  OSL^S  cm. 
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1.5.2.  Three  Degree  of  Freedom  Model 

The  SDOF  model  of  Figure  1.2  is  not  adequate  when  the  shafts  and  bearings  are 
compliant.  To  overcome  this  deficiency,  a  three  degree  of  freedom  (3-DOF) 
transverse -torsional  model  as  shown  in  Figure  1 .1  la  is  developed.  Equation  1.4  gives 
the  equations  of  motion  with  dimensionless  mass  [M],  damping  [C]  and  stiffness  [K] 
matrices  and  displacement  vector  { q)  given  as  follows: 


[M]  = 


1 

0 

0 


0 

nrigi  /  mci 
0 


0 

0 

m g2  /  mcj 


(1.14a) 


[C]  = 


Ch 

v/kh™cl 


1 

1 

-1 


1  -1 

(1  +  cbsl/ch)  _1 

“1  +  cbs2  '  ch  \ 


(1.14b) 


[K]  = 


1 

1 

-1 


1  -1 

(1  +  kbsl/kh) 

-1  (l  +  kbs2 /kb) 


(1.14c) 


{q}  =  [u,  ygI,  >’g2 ] 


(1.1 4d) 


8gl  ^g2®g2 
u  =  — 2 - £ — ; 

2  2dgl 


i  =  1,2. 


(1 . 14e,f) 


Here,  kbsj  and  kbs2  are  equivalent  lateral  stiffnesses  representing  shaft  and  bearing 
flexibilities,  and  cbsj  and  cbs2  are  equivalent  viscous  damping  values.  This  model  can 
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Figure  1.11.  Reduced  order  analyt 


accurately  predict  all  three  modes  of  interest  as  evident  from  Figure  2.8  and  Table  1.3 
where  its  predictions  are  compared  with  the  results  of  FEM. 

When  the  system  is  connected  to  the  motor  and  load  inertias,  then  the  in  between 
torsional  stiffnesses  Ki  and  K2  should  be  compliant  enough  to  be  abie  to  neglect  the 
effects  of  the  motor  and  load  inertias,  as  it  is  discussed  earlier  in  Section  1.4.4.  In 
summary,  the  3-DOF  model  shown  in  Figure  1.11a  can  be  used  to  describe  the 
dynamics  of  the  geared  rotor  system  when:  a)  the  shafts  and  bearings  are  compliant  and 
provided  the  shafts  are  short  such  that  higher  order  bending  modes  of  the  shafts  are  out 
of  the  frequency  range  considered,  and  b)  the  torsional  stiffnesses  of  the  connections  in 
between  the  motor  and  load  inertias  and  the  gear  box  are  sufficiently  compliant. 


1.5.3.  Six  Degree  of  Freedom  Model 

A  six  degree  of  freedom  (6-DOF)  model  as  shown  in  Figure  1.11b  can  be 
employed  to  represent  the  geared  rotor-bearing  system  when  the  effects  of  the  motor 
and  load  are  not  negligible  as  mentioned  in  Section  1.4.4.  Equations  of  motion  are  still 
given  by  equation  (1.4)  and  the  dimensionless  system  matrices  are  defined  as 


[Mi¬ 


n'd 


Ip/dgl 
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0 
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Igi/dgi 

0 

0 
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0 
0 

g2/dgi 
0 
0 
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0 
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iL/dfi 
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mgl 
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mg2 


(1.15a) 
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Table  1.3.  Comparison  of  typical  modes  obtained  by  FEM  and  3-DOF  models  for  set 
B  of  Table  1.1;  Lsj/dgj=0.3. 


Mode  \\i 

Displacement 

Vi 

FEM/3-DOF 

VU 

FEM/3-DOF 

Vm 

FEM/3-DOF 

-0.012/-0.01 1 

-0.001/0.0 

-0.210/-0.246 

0g2 

0.012/0.011 

-0.001/0.0 

0.210/0.246 

ygi 

1.0/ 1.0 

1. 0/1.0 

-1.0/- 1.0 

yg2 

-1.0/- 1.0 

1. 0/1.0 

1. 0/1.0 
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0  0  0 
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_£h 

2 

2 

C2 

Chvg 
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2 

L, 

Cj. 
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0 

(ch  +cb*l>  ”* Ch 

(Ch  +  cb»2>. 


(1.15b) 


[K]  = 
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i  +ih( 
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khv 


dgl 


hj_g 
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khv 


-symmetric  - 
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khvg 
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(kh  +  kb*l ) 
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~  2 
khvg 
2 


~kh 

(kh  +kbs2) 


(1.15c) 


{q}  —  ^9p,  ®gi»  ®g2>  ^l>  ygi’  y  ] 


(1.1 5d) 


For  data  set  B  of  Table  1.1  predictions  yielded  by  the  6- DOF  model  are  compared 
with  those  by  FEM  as  shown  in  Figure  2.9  and  Table  1.4.  Based  on  these  results,  tt 
can  be  concluded  that  the  6-DOF  model  is  accurate  enough  to  predict  the  natural  modes 
of  the  system.  Accordingly,  6-DOF  model  must  be  employed  when  the  effects  of  the 
motor  and  load  are  not  negligible. 
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Tabic  1.4.  Comparison  of  typical  modes  obtained  by  FEM  and  b-DOF  models  for  set 
B  of  Table  1.1;  L^/d^O-S,  Ip/Igl=5. 


Mode  y 

Displacement 

Vi 

FEM/6-DOF 

Vn 

FEM/6-DOF 

Vie 

FEM/6-DOF 

0P 

0.079/0.079 

0.002/0.0 

-0.018/-0.013 

-1.0/- 1.0 

-0.005/0.0 

1. 0/1.0 

0g2 

1. 0/1.0 

-0.005/0.0 

-1.0/- 1.0 

eL 

-0.079/-0.079 

-0.002/0.0 

0.018/0.013 

ygi 

0.006/0.006 

1. 0/1.0 

0.391/0.451 

yg2 

-0.006/-0.006 

1. 0/1.0 

-0.391/-0.451 
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1.6.  CONCLUSION 

In  this  chapter,  a  finite  element  model  to  investigate  the  dynamic  behavior  of 
linear  time-invariant  geared  rotor  systems  has  been  developed.  The  transverse 
vibration  of  the  system  associated  with  shaft  and  bearing  flexibilities  and  the  dynamic 
coupling  between  the  transverse  and  torsional  vibrations  due  to  gear  mesh  have  been 
considered.  Natural  modes  of  the  system  have  been  identified  and  forced  vibration 
response  due  to  both  low  frequency  external  and  high  frequency  internal  excitations 
have  been  determined.  Reduced  order  analytical  models  of  the  geared  rotor  bearing 
system  have  also  been  developed.  Three  different  linear  time-invariant  models  (SDOF, 
3-DOF  and  6-DOF)  have  been  suggested  to  represent  the  geared  system.  By 
comparing  results  with  FEM  predictions,  it  has  been  shown  that  such  reduced  order 
linear  models  are  reasonably  accurate.  Therefore  these  models  will  be  extended  in 
Chapters  II,  III  and  IV  in  analyzing  the  effects  of  system  non-linearities  and  time- 
varying  mesh  stiffness. 


36 


CHAPTER  II 


NON-LINEAR  DYNAMIC  ANALYSIS  OF 
A  SPUR  GEAR  PAIR 

2.1.  INTRODUCTION 

2.1.1.  Excitation  Types  and  Backlash 

The  focus  of  this  chapter  is  on  the  backlash  non-linearity  as  excited  primarily  by 
the  transmission  error  between  the  spur  gear  pair.  A  gear  pair  is  bound  to  have  some 
backlash  which  may  be  designed  to  provide  adequate  lubrication  and  eliminate 
interference  due  to  manufacturing  errors.  Backlash-induced  torsional  vibrations  may 
cause  tooth  separation  and  impacts  in  unloaded  or  lightly  loaded  geared  drives.  Such 
impacts  result  in  intense  vibration  and  noise  problems  and  large  dynamic  loads, 
which  may  affect  reliability  and  life  of  the  gear  drive  [16,17].  Excitation  mechanisms 
can  be  grouped  as  follows: 

A.  External  Excitations:  This  group  includes  excitations  due  to  rotating  mass 
unbalances,  geometric  eccentricities,  and  prime  mover  and/or  load  torque  fluctuations 
[18].  Although  mass  unbalances  and  geometric  eccentricities  can  be  reduced  through 
improved  design  and  manufacturing,  torque  fluctuations  are  not  easy  to  eliminate 
since  they  are  determined  by  the  characteristics  of  the  prime  mover  (piston  engines, 
dc  motors  etc.)  and  load  [19].  Such  excitations  are  typically  at  low  frequencies  f2T 

which  are  the  first  few  multiples  of  the  input  shaft  speed  f 2S .  Practical  examples 
include  rattle  problems  in  lightly  loaded  automotive  transmissions  and  machine  tools 
[19,20]. 
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B.  Internal  excitations;  This  group  includes  high  frequency  excitations 

caused  by  the  manufacturing  related  profile  and  spacing  errors,  and  the  elastic 
deformation  of  teeth,  shafts  and  bearings.  Under  the  static  conditions,  all  such 
mechanisms  can  be  combined  to  yield  an  overall  kinematic  error  function  known  as 
"the  static  transmission  error"  e(0  [17,18].  This  error  is  defined  as  the  difference 
between  the  actual  angular  position  of  the  driven  gear  and  where  it  would  be  if  the 
gears  were  perfectly  conjugate  [17,18,21-23].  In  gear  dynamic  models,  e(0  is 
modeled  as  a  periodic  displacement  excitation  at  the  mesh  point  along  the  line  of 
action  [15,24-26]  and  its  period  is  given  by  the  fundamental  meshing  frequency 
Qh=NQs  where  N  is  the  number  of  teeth  on  the  pinion  Practical  examples 

include  steady  state  noise  and  vibration  problems  in  automotive,  aerospace, 
industrial,  marine  and  appliance  geared  systems. 

2.1.2.  Literature  Review 

Experimental  studies  on  the  dynamic  behavior  of  a  spur  gear  pair  with  backlash 
started  almost  30  years  ago  and  still  continue  [27-29].  As  one  of  the  better  examples 
of  such  experiments,  Munro  [27]  developed  a  lightly  damped  (damping  ratio  £®0.02) 
four-square  test  rig  to  measure  the  dynamic  transmission  error  of  a  spur  gear  pair. 
He  used  high  precision  gears  with  rigid  shafts  and  bearings,  and  showed 
experimentally  that  the  tooth  separation  takes  place  when  the  mean  load  is  less  than 
the  design  load.  Dynamic  transmission  error  versus  speed  curves  were  ploned  to 
illustrate  the  steady  state  response  and  the  jump  phenomenon.  Kubo  [28]  measured 
the  dynamic  tooth  stresses  using  a  similar  set-up  in  order  to  calculate  the  dynamic 
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factors.  He  also  observed  a  jump  in  the  frequency  response  of  the  gear  pair  with 
backlash  even  though  the  test  set-up  was  heavily  damped  (£=0.1). 

Such  experimental  studies,  though  limited  in  scope,  have  clearly  shown  that  the 
gear  pair  dynamics  can  not  be  predicted  with  a  linear  model  -  see  Ozguven  and 
Houser  [14]  for  a  detailed  review  of  the  linear  gear  dynamic  models  as  available  in  the 
literature.  Although  most  of  the  non-linear  mathematical  models  used  to  describe  the 
dynamic  behavior  of  a  gear  pair  are  somewhat  similar  to  each  other,  they  differ  in 
terms  of  the  excitation  mechanisms  considered  and  the  solution  technique  used.  For 
instance,  a  large  number  of  studies  have  focused  on  the  rattle  problem  in  lightly 
loaded  geared  drives  which  are  excited  by  the  low  frequency  external  torque 
excitations  [30-35],  A  few  investigators  have  included  the  static  transmission  error 
excitation  in  the  non-linear  models  [24,35-37]. 

The  gear  backlash  non-linearity  is  essentially  a  discontinuous  and  non- 
differentiable  function  and  i;  represents  a  strong  non-linear  interaction  in  the 
governing  differential  equation.  This  issue  has  been  discussed  by  Comparin  and 
Singh  [33]  and  they  have  concluded  that  most  of  the  solution  techniques  available  in 
the  literature  can  not  be  directly  applied  to  examine  this  problem.  Most  of  the  gear 
dynamic  researchers  recognized  this  problem  implicitly,  and  therefore  employed 
either  digital  or  analog  simulation  techniques  [19,24,29,35-38].  For  instance, 
Umezawa  et.al  [29],  Yang  and  Lin  [32]  and  Ozguven  and  Houser  [24]  have  solved  a 
one  degree  of  freedom  torsional  model  of  the  gear  pair  using  numerical  techniques. 
Lin  et.al.  [38]  included  motor  and  load  inertias  in  a  three  degree  of  freedom  torsional 
model.  Kucukay  [35]  has  developed  an  eight  degrees  of  freedom  model  to  include 
the  rocking  and  axial  motions  of  the  rigid  shafts.  In  most  of  these  studies,  with  the 
exception  of  Umezawa's  analysis  [29]  which  did  not  include  any  backlash,  a 
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discontinuity  has  been  seen  in  the  frequency  response  characteristics.  But  many 
investigators  have  typically  joined  two  discrete  points  to  show  a  broad  jump  in  the 
frequency  response  curve  [24,35,38]  Some  of  these  problems  have  been  due  to  the 
numerical  simulation  techniques  which  may  not  work  or  may  result  in  misleading 
answers  if  not  employed  properly.  Such  difficulties  have  been  found  by  Comparin 
and  Singh  [33],  Singh  et  al.  [34]  and  Gear  [39-40]  but  are  yet  to  be  resolved  or 
addressed  by  the  gear  dynamics  researchers.  Accordingly,  one  of  the  major 
objectives  of  this  chapter  is  to  examine  whether  numerical  simulation  techniques  can, 
in  fact,  be  used  to  predict  the  dynamic  response  completely,  and  what  precautions  one 
must  take  to  develop  such  a  mathematical  model.  Since  Comparin  and  Singh  [33] 
and  Singh  et  al.  [34]  have  examined  the  external  excitation  problem,  in  this  chapter 
we  focus  mainly  on  the  internal  excitation  and  see  whether  the  numerical  simulation 
technique  can  be  made  to  work  for  the  prediction  of  the  non-linear  frequency 
response  characteristics. 

A  few  researchers  have  attempted  to  obtain  the  analytical  solutions  for  a  gear 
pair  problem,  based  on  the  piece-wise  linear  techniques  which  divided  the  non-linear 
regime  into  several  linear  regimes  [41-43].  For  instance,  Wang  [41,42]  has  used  two 
and  three  degree  of  freedom  torsional  models  with  backlash,  and  assumed  that  the 
gear  teeth  are  rigid  and  the  driven  gear  has  an  infinite  inertia.  The  governing 
equations  have  been  solved  using  the  piece-wise  linear  technique.  It  should  be  noted 
that  the  piece-wise  linear  technique  gives  only  solutions  for  the  equivalent  linear 
systems  and  one  typically  may  have  difficulties  in  combining  such  solutions  [43]. 
Comparin  and  Singh  [33]  overcame  these  problems  by  employing  the  harmonic 
balance  method  (HBM)  and  constructed  analytical  solutions  for  the  non-linear 
frequency  response  characteristics  of  a  gear  pair  with  backlash  as  excited  by  the 
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external  torque.  In  this  chapter,  we  will  use  the  same  technique  to  examine  the 
internal  excitation  problem,  and  compare  results  with  digital  simulation  and 
experimental  studies.  Further  literature  review  is  included  in  subsequent  sections. 


2.2.  PROBLEM  FORMULATION 
2.2.1.  Physical  Model 

A  two  degree  of  freedom  semi-definite  model  of  the  spur  gear  pair  with  rotary 
inertias  Ig,  and  Ig2  and  base  circle  diameters  dgl  and  dg2  as  shown  in  Figure  2.1  is 

considered  here.  The  shafts  and  bearings  are  assumed  to  be  rigid.  The  gear  mesh  is 
described  by  backlash  of  2b  and  by  a  :imc  invariant  mesh  stiffness  kh  *  kh(t) 
when  in  contact  and  viscous  damping  cfa.  The  equations  of  torsional  motion  of  the 

gear  pair  as  shown  in  Figure  2.1  are 


-T..<n 


(2.1a) 
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g2 


_2 

dt 


d  ,c.  fd  ,  d0  , 

g  2  h  gl _ g^ 

2  V  2  dF 


dg2 d9g2  del 

2  dF  dF  J 


(2.1b) 
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Figure  2.1.  Gear  pair  model 
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where  T  (t)*T  .  +  T  (t)  and  T  It )  =  T  +  T  (t )  are  torques 

gl  glm  gia  g  2  g2m  g2a 

on  pinion  and  gear  and  f  is  a  non-anaJytical  function  essentially  describing  the  mesh 

elastic  force  as  shown  in  Figure  2.1.  Here,  output  torque  fluctuation  T  _  (0  will 

K 

be  neglected  to  simplify  the  dynamic  problem,  i.e.  =  Tg2m.  Equations 

(2.1a)  and  (2.1b)  can  be  reduced  to  one  equation  in  terms  of  q(0  which  is  defined 
as  the  difference  between  the  dynamic  transmission  error  x(t)  and  the  static 
transmission  error  e(t ) . 

^  +  Ch  f  +  M«9)  “  +  F*T(i)  '  ^  W  :  (2'2a> 

q(r)=«D-?(0  =  ^fi-e1!l(r)--|ieg2(r)-e(r) ; 


m  ,  =  7 
cl  / 


u  a 

gl 

,  41  , 

V  gl 


I2' 

g2 


41 


g2  J 


41 

mc2  =  T 


li 

gl 


(2.2c, d) 


— 

m 


iV.  iV  .  p  (n  2m«T«u<n 
d8l  =  *,2  '  *T  "  "\2dgl 


(2.2e,f) 


f(q(D) 


f(q(0) 


q(D  -  b; 

0; 

q(D  +  b; 
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where  mcl  is  the  equivalent  gear  pair  mass,  Fm  is  the  average  force  transmitted 
through  the  gear  pair,  FaT(t )  is  the  fluctuating  force  related  to  the  external  input 

torque  excitation  and  f(q(T))  is  the  nonlinear  displacement  function.  Equation 
(2.2a)  is  nondimensionalized  by  letting  q(t")  =  q(t~)  /  b,  05 n  =  \/^h  ^  m  d  * 
t  =  (o„r  and  ^=ch/(2mcl(0nj.  Now,  consider  harmonic  excitation  for  both 

e(D  and  F^t")  as  e(T)  =  e  sin(flhr  +  <t>h),  FaT(0  =  F  aT  sin  (  ft  Tt  +4>T) 
where  fth  and  ftT  are  the  fundamental  excitation  frequencies  of  internal 

displacement  and  external  torque  fluctuations,  respectively.  Further,  define 
dimensionless  excitation  frequencies  fth  =  fth/o>n  and  ftT  =  ftT/(on, 

dimensionless  external  mean  load  Fm  =  Fm/bkh,  amplitudes  of  the  dimensionless 
internal  (F^t))  and  external  alternating  forces  (Fe(t))  FaT=F#T/bkh  and 
F  ^  =  e  /  b  and  nonlinear  displacement  function  f(q(t))  to  yield  the  following 
governing  equation  of  motion. 


q(t)  +  2^q(t)+ f(q(t))  =  F(t)  (2.3a) 

F(t)=Fm  +  FaTsin(QTt  +<|>T)+Fabn2h  sin(ftht  +  4>h)  (2.3b) 


f(q(t))  = 


f(q(0) 

b 


q( t )  —  1 ;  q(t)  >  1 

•  0;  -  1  <  q(t)  <  1 

q(0  +  1;  q(t)  <  -  1 


(2.3c) 
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2.2.2.  Scope  and  Objectives 

When  only  external  forces  excite  the  system,  i.e.  Fj(t)=0,  equation  (2.3a) 
reduces  to 

q(t)+  2£q(t)+  f ( q( t ) )  =  Fm  +  Fe(t)  =  I  fll  +  FftX  sin  +  <(>T)  (2.4) 

This  equation  has  beer,  solved  both  analytically  and  numerically  by  Comparin  and 
Singh  [33].  Conversely,  no  analytical  solution  is  available  when  the  system  is 
excited  by  internal  static  transmission  error  at  the  mesh  frequency  which  is 
considerably  higher  than  The  governing  equation  is  given  by  substituting 
Fe(t)=0  in  equation  (2.3a);  note  that  the  external  mean  load  Fm  is  not  equal  to  zero. 

q(t)  +  2Cq(0  +  f(q(t))  =  Fm  +  F.(t)  =  Fm  +  F^ft^  sin(Qht  +  <j>h)  (2.5) 

Both  equations  (2.4)  and  (2.5)  include  the  clearance  non-linearity.  While 
equation  (2.4)  represents  the  conventional  represc  Nation  of  the  vibro-impact  problem 
[33,34],  equation  (2.5)  is  more  applicable  to  the  clearance  problems  in  built-up 
assemblies  where  the  excitation  is  generated  by  the  kinematic  errors.  This  chapter 
focuses  on  the  steady  frequency  response  characteristics  of  equation  (2.5)  which 
represents  a  gear  pair  with  backlash  as  excited  harmonically  by  the  static  transmission 
error  excitation  e(F)  orF^.  Specific  objectives  of  this  chapter  are  as  follows: 

1.  Solve  equation  (2.5)  numerically  to  resolve  various  modeling  issues  such  as 
the  existence  of  multiple  solutions,  subharmonic  resonances  and  chaos. 
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2.  Construct  analytical  solutions  to  equation  (2.5)  using  the  harmonic  balance 
method  (HBM)  which  has  been  applied  successfully  to  solve  equation  (2.4)  by 
Comparin  and  Singh  [33]. 

3.  Compare  digital  simulation  and  harmonic  balance  techniques  and  establish 
the  premises  under  which  the  jump  phenomenon  can  be  predicted. 

4.  Perform  parametric  studies  in  order  to  understand  the  effects  of  FnvFah  and 

A 

4  on  the  frequency  response.  Vary  the  force  ratio  FsT^/F^  which  is  a  measure  of 
the  load  on  the  gear  pair  and  compare  the  dynamic  behavior  for  lightly  and  heavily 
loaded  gears. 

5.  Validate  analytical  and  numerical  solution  techniques  by  comparing  these 
with  previous  experimental  studies  [27,28]. 

6.  Compare  the  frequency  response  characteristics  of  equations  (2.4)  and  (2.5), 
and  also  examine  the  possibility  of  finding  overall  response  when  both  external  and 
internal  excitations  are  applied  simultaneously. 

7.  Consider  the  periodic  static  transmission  error  excitation  case,  i.e. 

k  2 

F(t)=  Fm+F.(t)=Fm+  X  (Jfth)  Fahisin(jnht  + V5  °nly  ,hC  firS* 

j=  1 

three  (k=3)  harmonics  are  included. 

2.3.  DIGITAL  SIMULATION 

Clearance  or  vibro-impact  problems  in  single  degree  of  freedom  systems  have 
been  examined  by  a  number  of  investigators  whose  formulations  arc  similar  to 
equation  (4)  -  see  Comparin  and  Singh  [33]  for  a  detailed  review.  Moreover,  Shaw 
and  Holmes  [43]  and  Moon  and  Shaw  [44]  have  considered  an  elastic  beam  with  one 
sided  amplitude  constraint  subject  to  a  periodic  displacement  excitation,  and  have 
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shown  experimentally  and  numerically  that  the  chaotic  and  subharmonic  resonance 
regimes  exist.  Whiston  [45-47]  has  investigated  the  non-linear  response  of  a 
mechanical  oscillator  preloaded  against  a  stop.  He  has  solved  the  system  equation  for 
harmonic  excitation  by  digital  simulation  and  studied  the  existence  and  stability  of  the 
subharmonic  and  chaotic  responses  and  the  effect  of  preload  on  chaos.  Similarly, 
Ueda  [48]  has  solved  the  Duffing's  equation,  q  +  2£q  +  q3  =  F  sin  t ,  numerically 

and  defined  the  regions  of  different  solutions  on  a  £  versus  F  map.  According  to 
him,  the  existence  of  harmonic,  subharmonic  and  chaotic  responses  depends  on 
values  of  £  and  F,  and  multiple  steady  state  solutions  typically  exist.  Thompson  and 
Stewart  [49]  have  reviewed  the  available  literature,  with  focus  on  the  Duffing's 
equation.  It  should  be  noted  that  equation  (2.5)  is  different  from  the  non-linear 
differential  equations  considered  by  the  above  mentioned  studies.  Therefore, 
equation  (2.5)  must  be  studied  in  depth  as  the  results  of  the  other  non-linear  equations 
may  not  be  directly  applicable  to  our  case. 

First,  we  solve  the  governing  non-linear  differential  equation  (2.5)  numerically 
using  a  5th-6th  order,  variable  step  Runge-Kutta  numerical  integration  routine 
(DVERK  of  IMSL  [50])  which  is  suitable  for  a  strongly  non-linear  equation 
[33,39,40],  Second,  we  investigate  the  existence  of  chaos  and  subharmonic 
resonances.  Since  the  steady  state  response  of  the  system  due  to  the  sinusoidal 
excitation  is  of  major  interest,  it  is  necessary  to  run  the  numerical  program  for  a 
sufficient  length  of  time.  The  number  of  cycles  of  the  forcing  function  required  to 
reach  the  steady  state  depends  on  £. 

A  lightly  loaded  system  (F=Fm./Fah=0.5)  with  low  damping  (£=0.02)  is 
considered  as  the  first  example  case.  For  Fm=0.1  and  Fah=°-2.  the  response  q(t)  is 
computed  over  the  frequency  range  0  <  <  1.5.  Figure  2.2  shows  phase  plane 
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plots  q(t)  vs.  q(t)  for  different  Qh  values.  For  Qh=0.3,  all  transients  converge  to 
one  periodic  solution  at  the  fundamental  frequency  £ih  of  the  forcing  function 
irrespective  of  the  initial  conditions  q(0)  and  q(0) .  Therefore,  it  is  called  a  "period- 
one,  tp  "  attractor  where  tp=27t/Qh.  But,  in  the  case  of  Qh=0.5,  three  coexisting 
period-one  attractors  have  been  found  as  shown  in  Figure  2.2.  Here,  q(0)  and  q(0) 
define  three  steady  state  limit  cycle  solutions.  For  all  initial  conditions  given  by 
—  2  <  q(0 )  <  2  and  -  2  <  q(0 )  <  1 ,  a  map  of  the  domains  of  attraction  for  each 
steady  state  solution  is  obtained  in  the  map  in  Figure  2.3.  If  a  smaller  increment  is 
used  for  the  initial  conditions,  a  finer  resolution  will  be  obtained.  Hence  each  phase 
plane  plot  shown  in  Figure  2.2  is  strictly  governed  by  a  subset  of  the  initial 
conditions.  Similarly,  two  period-one  attractors  are  found  at  Qh=0.6.  At  Qh=0.7, 
besides  two  period-one  tp  attractors,  two  more  solutions  of  period  3tp  exist,  i.e. 
period-three  attractors.  But  only  a  period-two  2tp  attractor  is  seen  when  fth=0.8. 
Within  the  range  1.  0  <*  ft  h  <,  1.  5 ,  non-periodic,  steady  state  or  chaotic  response  is 

observed.  Figure  2.4  shows  the  chaotic  time  history  and  the  Poincare  map  (strange 
attractor)  at  ft^l.O.  These  results  are  qualitatively,  but  not  quantitatively,  similar  to 

the  studies  reported  on  the  Duffing's  equation  [48,49]  and  clearance  non-linearity 
(equation  (2.4)  type)  [33]. 

It  is  concluded  that  the  subharmonic  response  of  period  ntp  provided  n  *  1  and 
the  chaotic  response  ( tp  —►«>)  are  seen  in  the  gear  pair  only  for  a  certain  set  of 

parameters  Fm,  £ ,  £  and  fth.  It  must  also  be  noted  that  the  multi-solution  regions 
are  strongly  dependent  on  the  choice  of  q(0)  and  q(0).  Only  one  steady  state 
solution  can  be  found  via  digital  simulation  when  orJy  one  set  cf  initial  conditions  is 
chosen  at  a  given  fth,  and  the  rest  of  the  steady  state  solutions  are  not  predicted.  This 
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initial  displacement 


Domains  of  attraction  for  three  attractors;  *  no  impact  (Figure  2.2bl),  • 
single  sided  impact  (Figure  2.2b3),  °  double  sided  impact  (Figure  2.2b2); 

F  =0.5,  Fm=0.1,  ^=0.02  and  Qh=0.5. 
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results  in  an  incomplete  frequency  response  description.  These  issues  will  be 
discussed  further  in  Section  2.5. 

A 

A  heavily  loaded  system  with  F=2,  Fm=0. 1  and  the  same  amount  of  damping 

£=0.02  is  considered  as  the  second  example  case.  Figure  2.5  shows  the  phase  plane 
plots  for  the  same  values  of  which  are  ured  in  the  first  example  case.  Unlike  the 

first  example  case,  no  chaotic  responses  are  found  here.  All  of  the  solutions  are 
period-one  tp  type  for  all  Qb<1.0  and  Figure  2.6  shows  the  typical  domains  of 
attraction  at  £2b=0.6.  However,  within  the  range  1.  0  £  Q  b  S»  1.  5 ,  one  period-two 

A 

attractor  co-exists  with  the  period-one  orbit.  Hence  the  force  ratio  F  determines  the 
existence  of  the  chaotic  and  subharmonic  responses.  To  illustrate  this  point,  consider 
chaos  shown  in  Figure  2.4.  As  £  and  £  are  increased,  significant  changes  in  the 
response  are  observed  in  Figures  2.7  and  2.8.  A  transition  from  the  chaos  to  a 
period-two,  and  then  to  a  period-one  steady  state  solution  is  seen  when  £  is 
increased  from  0.75  to  1.0  and  then  to  1.5.  Similarly,  an  increase  in  £  to  0.05 
reduces  chaos  to  a  period-eight  attractor,  which  then  bifurcates  to  a  period-two  orbit 
at  £=0.1  as  shown  in  Figure  2.8.  Since  most  real  geared  systems  are  heavily  loaded 

A 

with  a  high  F ,  chaotic  and  subharmonic  responses  should  not  be  seen  under  the 
normal  driving  conditions.  This  issue  will  be  discussed  again  in  section  2.6. 

2.4.  ANALYTICAL  SOLUTION 

An  approximate  solution  for  equation  (2.5)  is  constructed  using  the  harmonic 
balance  method  (HBM).  Assume  that  q(t)  =  qm  +qa  sin(Obt  +  <|>r)  where  q^  and  qa 
are  mean  and  alternating  components  of  the  steady  state  response,  and  $  r  is  the 
phase  angle.  Here,  higher  harmonics  of  the  response  are  not  included  in  the  analysis. 
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Domains  of  attraction  for  two  attractors;  *  no  impact  (Figure  2.5c2),  •  single 

A 

sided  impact  (Figure  2.5c  1);  F  =2,  Fm=0.1,  £=0.02  and  ty,=0.6. 
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Figure  2.7.  Time  histories  of  q(t)  for  Fah=0.2,  £=0.02  and  a) 
c)F  =  l  .5, 
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Figure  2.8. 
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Time  histories  of  q(t)  for  Fm=0.1,  F  =0.5  and  a)  £=0.05,  b)  £=0. 
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The  quasilinear  approximation  to  the  nonlinear  function  f(q)  with  the  excitation 
F(t)  =  Fm  +  F.(t) =  Fm  +  F^O ^  sin  (£2  fat  +  is  in  the  form  [33,51] 


f(<l)  =Nm<lin+  +<*>r)  +  N»q®C0S(aht  (26) 


where  the  describing  functions  Nm,  Na  and  N*  are  defined  as: 


N 


m 


_ 1 

2^qm 


2  n 


/f(qm  +q«  sin  <p)d<p  ; 
o 


(2.7a) 


Na  = 


req, 


2  n 

lf(  q. 

0 


+  qa  sin  (p)sin  <pd<p  ; 


(2.7b) 


N 


* 

a 


i 

icq. 


2  n 

jf(qm  +qa  sin  <p)cos<pdq>  ; 
o 


<P=  *V+  4>r 


(2.7c, d) 


Equation  (2.3b)  is  substituted  into  equations  (2.7a),  (2.7b)  and  (2.7c)  to  obtain 


Nm=>-  2^ 


m 


[g(Y+)-g(Y-)l  ;  Na=  1-  i[h(Y+)+h(Y-)]  ;  (2.8a, b) 


*±  = 


l±q 


m 


(2.8c,d) 


where 
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g(Y)  = 


|(Ysin  *Y  +  V1  -Y2)  ;  lYl  S  1  . 

Iy|;  IyI >  i  ’ 


n(sin”1Y  +  Y*yi-Y2)  ;  WSI 
h(Y)  =  ]  -  1  ;  y<-1 

11;  Y>1 


(2.8e) 


(2.8f) 


Comparin  and  Singh  [33]  used  truncated  series  expansions  for  functions  g(Y)  and 
h(Y)  given  by  equations  (2.8e)  and  (2.8f).  They  stated  that  the  error  involved  in 
using  truncated  series  is  within  6  percent  when  only  the  first  two  terms  are  considered 
with  the  coefficient  of  the  second  term  adjusted  to  yield  the  actual  value  for  y  =  1  • 
Using  the  same  approach,  one  gets 

8(Y)£|(1+(^)y2);  h(Y)S^Y“('iTiLy);  M*1  (2.9a, b) 


and  obtains  the  following  frequency  response  by  substituting  equations  (2.6)  into 
equation  (2.5)  and  equating  the  coefficients  of  like  harmonics: 


qa  = 


F  J1! 

ah  b 


7(N»-ft2j2+ 


\  Q  m 


N. 


(2.10a,b) 


4>r  =  4>h- tan 


-  1 


'  2;nh  ' 


vNa-^'hy 


(2.100 
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Depending  on  the  damping  ratio  £  and  the  parameters  Fm-  Fah  and  (2  which 
define  the  excitation,  there  are  three  cases  at  which  different  solutions  are  obtained: 
(a)  no  impact  (no  tooth  separation),  (b)  single  sided  impacts  (tooth  separation,  but  no 
back  collision),  and  (c)  two  sided  impacts  (back  collision). 

Case  I:  No  Impact  :  The  tooth  separation  (impact)  is  not  observed  in  a  geared 
system  if  the  displacement  q(t)  lies  in  the  region  q(t)>l  all  the  time.  This  condition, 
shown  in  Figure  2.9  as  case  I,  can  be  described  mathematically  as 

|qm+qal>i  and  |qm  - qa| > i  (2.1  n 

Then,  for  no  impact  region  in  which  the  conditions  defined  by  equation  (2.1 1)  are 
satisfied,  the  describing  functions  are  given  by 

Nm  =  1  “  q~~  ;  N,  -  1  •  (2.12a,b) 


Substitution  of  equations  (2.12)  into  equations  (2.9)  yields  the  following  governing 
equations  for  no  impact  case 


2 

f  ci 

ah  h 


qal  I  Z~T:  ; 

J  0-«b)  +(2^h>2 


q  ,  =  F_  +  l 

Mml  m 


(2.13a, b) 


^hl  =<t>p-*an-  1 


i  -  n 


(2.13c) 
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Figure  -2.9.  Illustration  of  different  impact  regimes. 
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Hence,  whe.:  there  is  no  tooth  separation,  the  system  is  linear  and  the  mean  and 
alternating  components  of  the  response,  and  q^  are  uncoupled  from  each  other. 

The  transition  frequencies  from  the  no  impact  (linear)  regime  to  the  other  non-linear 
regimes  where  tooth  separation  occurs  are  found  by  substituting  q ,  —  |q  m  -  l|  into 
the  equations  (2.13a)  and  (2.13c)  as 


(2.14) 


where  £2(J  and  Ql2  are  the  transition  frequencies  from  no  impact  to  single  sided 
impact  regimes  below  and  above  resonance,  respectively. 

Case  II:  Single  Sided  Impact  ;  Mathematically,  the  single  sided  impacts  (tooth 
separation  without  back  collision)  are  observed  if 

qm  +  qa>1  sad  |q m  q al < 1  ■  tf.is) 

As  illustrated  in  Figure  2.9  as  case  II,  the  solution  remains  in  the  region  q(t)  >  -  1 . 
The  describing  functions  satisfying  equation  (2.15)  are  in  the  form 

Nm=l-  N„=l-i[l+h(Y-)]  (2.16a.b) 

*  rnTT 
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Substitution  of  equations  (2.10a)  and  (2.10b)  into  equations  (2.16a)  and  (2.16b) 
gives  the  describing  functions  for  the  single  sided  impact  case: 

Nm =  i+  2Ta-{i[l  +  (ii1)<Y-)2]-(Y+>}  ;  <217a) 

M  mil 

1  2(l~qmll)r,  r  4  -  Tt  ,21 

Na=2 - nq'^-  )(Y->  J  (2.17b) 

Then  equations  (2.9a),  (2.9b)  and  (2.9c)  with  the  describing  functions  given  by 
equations  (2.17a)  and  (2.17b)  define  the  response  of  the  system  in  the  single  sided 
impact  region.  In  this  case,  it  is  hard  to  find  closed  form  expressions  for  transition 
frequencies  from  single  sided  to  double  sided  impact  regions.  The  validity  of  the 
solutions,  obtained  by  solving  equations  (2.9)  and  (2.17)  iteratively,  should  be 
checked.  If  the  solution  does  not  satisfy  equation  (2.15),  then  single  sided  impacts 
are  not  seen  at  that  particular  frequency. 

Case  III:  Double  Sided  Impact  :  The  double  sided  impact  case  exists  if  q^,,  and 
qg  are  such  that  the  following  conditions  are  satisfied: 

qaHl-qm|  arid  qa>jl+qj  (2.18) 

Figure  2.9  (case  III)  shows  the  double  sided  impact  case  at  which  qa  is  large  enough 
when  compared  to  q^  so  that  back  collision  is  observed.  The  describing  functions 
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for  this  case  are  obtained  from  equations  (2.8a)  and  (2.8b)  under  the  condition  given 
by  equations  (2.19): 


Nm=  1- 


2(  re  —  2) 


rtq 


am 


1  ( 

[1+3qHl 

v  4  J| 

<  qam  /. 

(2.19a,b) 


The  solution  for  the  double  sided  impact  case  is  obtained  by  solving  equations 
(2.9)  and  (2.19)  numerically.  The  validity  of  the  solution,  again,  should  be  examined 
by  using  the  conditions  defined  by  equation  (2.18),  as  it  is  done  in  the  single  sided 
impact  case. 

The  solutions  of  all  three  regimes  are  combined  to  obtain  the  overall  frequency 
response  of  the  gear  pair.  Figure  2.10  illustrates  typical  q*  versus  Qh  and  q,^  versus 
flh  plots.  All  three  impact  regimes  are  shown  on  these  plots.  Also  we  note  the 
frequency  region  where  multiple  steady  state  iolutions  are  obtained. 

2.5.  COMPARISON  OF  TWO  SOLUTION  METHODS 

First,  we  validate  the  approximate  analytical  solutions  of  equation  (2.5)  by 
comparing  predictions  with  the  results  obtained  by  digital  simulation.  Again,  two 
example  cases  (a  heavily  loaded  system  and  a  lightly  loaded  system)  are  considered 
and  the  frequency  response  curves  versus  Ob  and  qm  versus  f2h  are  generated. 

A 

Figure  2.1 1  shows  the  frequency  response  for  the  heavily  loaded  system  with  F=2, 
Fm=0.1  and  £=0.05.  Numerical  and  analytical  results  agree  very  well  as  both  predict 

amplitudes  and  the  transition  frequencies  for  case  I  (no  impact)  and  case  II  (single 
sided  Impact)  regimes,  and  both  show  that  the  double  sided  impact  solutions  do  not 
exist.  These  results  demonstrate  that  the  analytical  solutions  are  indeed  correct  fcr  the 


63 


frequency 


Figure  2.10.  Typical  frequency  response  plots  for  a  gear  pair;  a)  qa  versus  £2h,  b)  q 
versus  Qh. 
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heavily  loaded  system.  Now  consider  the  lightly  loaded  system  as  the  second 
example  case  with  F=0.5,  Fm=0.1  and  £=0.05.  From  Figure  2.12  we  note  that 

while  case  I  and  case  II  solutions  as  yielded  by  both  solution  techniques  are  very 
close  to  each  other,  case  HI  regime  is  predicted  only  by  the  analytical  expressions. 
Why  is  case  ID  not  predicted  by  the  digital  simulation?  To  answer  this  consider  the 
following: 

a)  In  digital  simulation,  several  sets  of  q(0)  and  q(0)  must  be  tried  to  find  all 

of  the  steady  state  solutions  within  a  multi-valued  region,  for  instance 

0.6  <  ft  £  0.  7  in  Figure  2.11  and  0.4  £  ft.  £  0.  7  in  Figure  2.12.  When  q(0) 
h  n 

and  q(0)  are  kept  constant  for  each  ftj,  or  if  the  steady  state  solution  of  the  previous 
frequency  is  used  as  the  basis  for  the  initial  conditions  for  the  next  frequency 
considered,  only  one  of  the  solutions  can  be  found  while  missing  the  other(s). 
Figure  2.13  illustrates  this  point  as  the  steady  state  solution  of  the  previous  frequency 
is  used  here  for  the  initial  conditions  at  each  ftb.  Using  this  procedure,  case  II 

solutions  within  the  multi-valued  regions  are  missed.  Therefore,  whenever  digital 
simulation  is  used  to  solve  the  equation  of  motion,  dependence  of  the  steady  state 
solutions  on  q(0)  and  q(0)  must  be  taken  into  account  to  avoid  the  risk  of  obtaining 
an  incomplete  frequency  response  description. 

b)  In  the  case  of  Figure  2.12,  no  initial  conditions  governed  by  case  III  are 
found  by  the  digital  simulation  within  the  range  of  initial  conditions  -  2  <  q(0  )  <  2 
and  -2<q(0)<  2  we  have  considered.  Conversely,  the  analytical  method 
predicts  this  regime  as  the  issue  of  initial  conditions  is  irrelevant  here.  To  illustrate 
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this,  first  we  consider  the  case  of  Figure  2.3  with  F=0.5,  Fm=0.1,  £=0.02  and 
fth=0.5.  Figure  2.14a  shows  the  domain  of  attraction  governed  by  case  III  within 
the  range  -  2  <  q(0 )  <  2  and  -  2  <  q(0)  <  2  with  a  q(0)  or  q(0)  increment  of 

0.2.  Here,  almost  half  of  the  initial  conditions  considered  define  case  HI  solutions. 
Therefore  the  chance  of  having  a  double  sided  impact  solution  at  Oh=0.5  is 

considerably  high.  Now,  increase  the  damping  to  £=0.03.  Once  again  a  case  III 
regime  is  predicted  analytically,  but  the  number  of  initial  conditions  corresponding  to 
case  m  is  not  sufficiently  high  like  the  case  of  £=0.02.  In  the  digital  simulation,  the 
transient  solution  will  converge  to  case  I  or  case  II  solution  unless  the  initial  condition 
corresponds  to  one  of  those  shown  in  Figure  2.14b  for  case  IQ.  And  the  limiting 
case  is  reached  when  £=0.04,  while  holding  the  other  parameters  the  same;  now  no 
initial  conditions  corresponding  to  case  III  are  found.  Obviously  there  might  be 
initial  conditions,  out  of  the  range  considered,  which  correspond  to  the  analytical 
results.  Therefore  the  existence  of  case  III  solutions  should  be  checked  numerically 
by  searching  the  entire  range  of  initial  conditions  as  defined  by  the  physical 
considerations  of  the  system. 

2.6.  PARAMETRIC  STUDIES 

Frequency  response  amplitudes,  transition  frequencies  and  the  existence  of 
various  impact  regimes  depend  on  Fm,  F  and  £.  Therefore,  a  set  of  parametric 
studies  using  analytical  solutions  of  section  2.4  will  be  presented  here;  same  results 
can  be  duplicated  by  the  digital  simulation  technique.  First  we  examine  the  effect  of 

A 

Fm  and  F  wliile  holding  the  damping  ratio  £  equal  to  0.05.  Figures  2.15  and  2.16 

A 

compare  results  for  four  different  F  values  obtained  by  varying  Fah  for  a  given 

A 

Fm=0.1.  Both  gears  maintain  complete  contact  with  each  other  when  F  is  very 
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Comparison  of  theory  with  numerical  simulation;  Fm=0.1, 
a)qa  versus  £2h,  b)  q^  versus  l2h. 


.5 


=2.  £=0.05; 


6 ' 


amplitude  amplitude 


(a) 


Figure  2.12. 


frequency 


(b) 


frequency 


Comparison  of  theory  with  numerical  simulation;  Fm=0.1 ,  F  =0.5,  ^=0.05, 
a)  qg  versus  Qh,  b)  qm  versus  ftb. 


68 


amplitude  amplitude 


(a) 


frequency 


Figure  2.13.  Comparison  of  theory  with  numerical  simulation  when  initial  conditions  are 
kept  constant;  Fm=0.1,  £=0.05;  a)  F=2,  b)  F=0.5. 
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large,  say  F=10  in  Figures  2.15a  and  2.16a.  Consequently  the  dynamic  system  is 
linear  and  the  mean  q^  and  alternating  qa  components  of  the  torsional  motion  are 

A 

uncoupled  as  expected.  However,  for  F  =2,  no-impact  regime  (case  I)  can  not  cover 
the  whole  frequency  range  and  a  region  around  the  resonant  frequency  is  dictated  by 
the  single  sided  impacts  (case  II).  Away  from  the  resonance,  over  the  range  given  by 
f2jj>0.8  and  f2j,<0.7,  solutions  are  single-valued,  whether  case  I  or  case  II  solutions. 
Conversely,  near  the  resonance  given  by  0.7<f2h<0.8  dual-valued  solutions  are  seen 
in  Figures  2.15b  and  2.16b.  And  sudden  discontinuities  are  observed,  i.e.  a  jump 

A 

up  at  nh=0.8  and  a  jump  down  at  fth=0.7.  When  F  is  reduced  further,  double 
sided  impact  (case  III)  solutions  start  to  appear  in  Figures  2.15c  and  2.16c.  Now,  the 
frequency  range  of  interest  can  be  divided  into  the  following  four  regions:  i)  case  I, 
Cth<0.4;  ii)  case  I  and  case  III ,  0.4<f2h<0.6;  iii)  all  three  cases,  0.6<flh<0.625;  and 

A 

iv)  only  case  II,  f2h>0.625.  Region  (ii)  disappears  at  F>0.5. 

A 

Next,  we  vary  F  by  changing  Fm  for  a  given  Fah  =0  .05.  This  set  of 
parameters  does  not  yield  any  double  sided  impacts  (Case  ID).  Again,  the  tooth  pair 
does  not  lose  contact  when  F  is  large  enough,  say  ^£10  and  only  the  linear 

A  A 

solutions  exist  as  shown  in  Figure  2.17a.  But  when  F  is  lowered  to  F  =4,  2  and  1, 
the  response  is  non-linear  which  is  composed  of  cases  I  and  II  as  shown  in  Figures 
2.17b-d. 

Next  we  examine  the  role  of  damping  ratio  £  on  the  frequency  response  in  Figure 
2.18,  given  Fm=0.1  and  Fah=2.  Double-sided  impacts  are  found  at  a  low  damping 

value  £=0.025  as  shown  in  Figure  2.18a.  When  £  is  increased  to  0  05,  case  III 
solutions  no  longer  exist  and  case  I  and  case  II  solutions  define  the  frequency 
response  completely.  The  jump-up  and  jump-down  transition  frequencies  in  Figure 
2.18b  are  distinctly  apart.  With  an  increase  in  £,  transition  frequencies  approach  each 
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Figure  2.15.  Frequency  response  qa  for  different  alternating  loads;  Fm=0.1,  £=0.05, 
a)  F  =10,  b)  F  =2,  c)  F  =  l,d)  F=0.5. 
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Figure  2.17.  Frequency  response  qa  for  different  mean  loads;  Fah=0.05,  £=0.05,  a) 

A  AAA 

F  =10,  b)  F  =4,  c)  F  =2,  d )  F  =  1 . 
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other  which  narrows  the  dual-valued  solution  region.  At  £=0.1,  the  frequency 
response  is  single  valued  and  the  jump-up  and  jump-down  transitions  in  Figure  2. 1 8d 
take  place  almost  at  the  same  frequency.  Finally,  like  the  linear  systems,  the  damping 
ratio  also  lowers  the  amplitudes  in  the  non-linear  resonance  regime. 

2.7.  EXPERIMENTAL  VALIDATION 

Our  analytical  solutions  of  Section  2.4  will  be  compared  with  the  experimental 
results  of  Munro  [27]  who  used  a  four-square  test  rig  to  measure  the  dynamic 
transmission  error  x(t)of  a  spur  gear  pair.  High  precision  spur  gears  with 
manufacturing  errors  much  smaller  than  tooth  deflections  were  selected.  Pinion  and 
gear  were  identical  with  32  teeth,  face  width  of  12.7  mm  and  diametral  pitch  of  4. 
Tooth  profile  modifications  were  applied  to  obtain  a  minimum  (but  not  zero)  e(t )  at 
the  design  load  of  3780  N.  Other  components  of  the  set-up  including  shafts, 
bearings  and  casing  were  made  as  rigid  as  possible  in  order  to  simulate  the 
configuration  shown  in  Figure  2.1.  x(t )  was  measured  for  a  range  of  gear  mesh 
excitation  frequencies  under  different  mean  loads  Fm  .  Some  of  the  key  parameters 

were  not  specified  by  Munro  127],  For  example,  it  was  stated  that  some  additional 
inertias  were  added  to  the  gears  to  shift  the  primary  resonant  frequencj  within  the 
operational  speed  range,  but  the  specific  values  of  such  inertit.r,  were  not  given.  It 
was  also  reported  that  the  damping  ratio  £  varies  with  load  "in  a  random  manner". 

Also,  backlash  was  not  measured  or  reported.  Therefore,  in  our  study,  we  estimate 
the  damping  ratio  (£=0.017.5)  and  the  resonant  frequency  o>  t)  by  considering  the 

design  load  case  a :  which  onlv  linear  tmavior  is  seen;  time-invariant  mesh  stiffness  is 
assumed  in  the  model.  The  same  vaiue  of  C  is  used  at  each  discrete  load  F_,  and  a 

backlash  value  of  2b=0. 1  mm  is  assumed  in  our  model 


7'» 


Figures  2.19  to  2.22  compare  measured  and  predicted  frequency  response  of 
the  dynamic  transmission  error  x(t ) ,  on  a  peak  to  peak  basis,  at  the  design  load  and 
at  1/4,  1/2  and  3/4  of  the  design  load,  respectively.  As  shown  in  Figure  2.19,  both 
analytical  and  experimental  results  indicate  that  teeth  maintain  contact  at  the  design 
load  and  hence  the  system  behaves  as  a  linear  system  in  spite  of  the  backlash.  This  is 
because  the  static  transmission  error  e(t  )  is  minimum  at  the  design  load  which 
results  in  a  large  force  ratio  Fm/Fa,  say  50.  However,  when  the  mean  load  Fm  is 

lowered  to  3/4  of  the  design  load  which  corresponds  to  a  larger  static  transmission 
error  e(t  ),  tooth  separation  takes  place.  Consequently  we  note  a  jump  in  the 
frequency  response  a  ;  shown  in  Figure  2.20  for  both  analysis  and  experiment.  This 
jump  phenomenon  is  more  noticeable  at  1/2  and  1/4  of  the  design  loads  as  shown  in 
Figures  2.21  and  2.22,  respectively.  Our  theory  matches  very  well  with  Munro's 
experimental  results  [27]  in  spite  of  the  lack  of  knowledge  of  some  relevant 
parameters. 

As  the  second  example  case,  the  experimental  results  of  Kubo  [28]  are 
considered  and  compared  with  our  theory.  Experimental  results  and  the  relevant 
system  parameters  are  extracted  from  a  recent  paper  by  Ozguven  and  Houser  [24]  and 
the  excitation  e(t )  is  calculated  using  a  spur  gear  elastic  model  [21].  Kubo 
designed  a  four-square  spur  gear  test  rig  which  was  heavily  damped  (£=0.1).  He 
measured  dynamic  root  stresses  od  and  then  estimated  the  dynamic  factor  Ds  as 
Ds=cid/crs  where  os  is  the  static  tooth  root  stress.  However,  in  several  other  studies 
[16]  the  dynamic  factor  is  defined  as  the  ratio  of  the  dynamic  mesh  load  F^  to  static 
mesh  load  Fm,  given  by  D L  =  F d  /  F  m  =  F /  F m  =  (2£q(t )  +  f(q(t ) )  /  Fm 
Note  that  DL  is  equal  to  Dg  when  the  change  of  the  moment  arm  due  to  a  change  in 
contact  point  is  neglected.  Figure  2.23  shows  the  envelope  of  dynamic  factor  Ds 
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Figure  2.21.  Comparison  of  theory  with  Munro's  [27]  experiments  at  1/2  of  design 
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Figure  2.22.  Comparison  of  theory  with  Munro's  [27]  experiments  at  1/4  of  desig 
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Comparison  of  theory  with  experimental  results  of  Kubo  (extracted  from 
reference  [24]). 
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measured  for  eight  different  teeth  pairs  and  DL  spectrum  computed  using  the 
analytical  solutions  of  Section  2.4.  Here,  at  most  of  the  frequencies,  our  predictions 
are  not  within  the  meastued  envelope,  but  the  transition  frequency  and  the  amplitude 
at  the  jump  discontinuity  are  predicted  accurately.  There  might  be  several  reasons  for 
the  discrepancy  including  the  usage  of  Dj_  instead  of  Ds,  the  validity  of  the 

computer  model  used  to  calculate  e(f),  insufficient  knowledge  of  some  system 
parameters  such  as  profile  modifications,  and  the  assumptions  made  in  developing 
our  theory  such  as  the  time  invariant  mesh  stiffness.  Therefore  the  experimental  data 
is  not  exactly  analogous  to  our  analytical  predictions  even  though  a  satisfactory 
agreement  is  found  Obviously  our  theory  should  be  refined  in  order  to  obtain  a 
better  agreement  with  Kubo’s  data;  for  instance,  time-varying  mesh  stiffness  may  be 
included. 

2.8.  COMPARISON  OF  EXCITATIONS 
2.8.1.  Internal  versus  External  Excitation 

First,  we  assume  that  only  one  type  of  excitation  exists  at  a  time  and  compare 
the  frequency  response  characteristics  of  the  system  due  to  the  internal  static 
transmission  error  sinusiodal  excitation  F^t)  given  by  equation  (2.5)  with  the  external 
sinusoidal  excitation  Fe(t)  given  by  equation  (2.4).  The  comparison  is  based  on  the 
analytical  solutions  which  are  constructed  in  Section  2.4  for  equation  (2.5)  and  in 

Reference  [33]  for  equation  (2.4).  In  the  case  of  internal  excitation,  the  amplitude  of 

2 

the  alternating  force  has  a  £2  h  term  which  makes  the  alternating  force  amplitude 

2 

frequency  dependent.  This  amplitude  is  smaller  than  Fah  forf2n<1  and 

greater  than  Fah  when  Hence  the  overall  alternating  force  amplitude  ratio 
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2  A 

Fm  /  varies  with  ftb  even  though  F  =  Fm  /  F^  is  kept  constant.  In  the 

A 

case  of  external  excitation,  the  force  amplitude  ratio  F=Fm/FaT  is  frequency 
independent. 

A 

Figures  2.24  and  2.25  compare  the  frequency  responses  for  four  values  of  F 
given  Fm=0.l,  ftb=ftr=ft  and  £=0.05.  When  F  is  sufficiently  large,  say  F  £10, 

only  case  I  solutions  corresponding  to  the  no-impact  (case  I)  exist.  Accordingly  both 
excitations  result  in  the  linear  system  responses  which  are  close  to  each  other  for 
ft<l;  but  for  ft>l  there  is  a  difference  which  grows  with  increasing  frequency  ft  as 

A 

shown  in  Figure  2.24a.  When  F  is  reduced  to  2,  both  case  I  and  case  II  regimes 
exist;  the  difference  between  two  excitations  is  again  significant  at  higher  ft  as  shown 
in  Figure  2.24b.  The  transition  frequencies  also  differ,  and  a  larger  range  of  dual¬ 
valued  solutions  is  seen  for  equation  (2.5).  Case  III  solutions  are  witnessed  at  lower 

A 

values  of  F  (1  and  0.5)  as  shown  in  Figures  2.25a  and  2.25b.  While  equation  (2.5) 
always  has  a  case  I  regime  at  low  ft,  equation  (2.4)  does  not  produce  case  I  at 

A 

F  =0.5.  Another  important  observation  from  Figures  2.24  and  2.25  is  that  up  to  two 
steady  state  solutions  are  seen  in  the  external  torque  excitation  case;  conversely  as 
many  as  three  solutions  are  found  for  the  geared  system  excited  by  the  static 
transmission  error. 

2.8.2.  Periodic  or  Combined  Excitation 

The  approximate  analytical  solutions  given  in  Section  2.4  are  constructed  only 
for  a  single  harmonic  internal  excitation  term  e(0  or  Fj(t).  However,  in  the  real 
geared  systems,  e(t  )  or  Fj(t)  is  a  periodic  which  can  be  represented  by  a  fourier 
series  of  the  fundamental  frequency  fth.  Therefore,  it  is  necessary  to  consider  the 
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Figure  2.24. 


Comparison  of  the  frequency  response  due  to  the  external  torque  excitation 
given  by  Comparin  and  Singh  [33]  and  the  static  transmission  enor  excitation, 
£=0.05,  Fm=0.l,Qh=QTand  a)  F=i0,b)  F=2. 
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Figure  2.25.  Comparison  of  the  frequency  response  due  to  the  external  torque  excitation 
given  by  Comparin  and  Singh  [33]  and  the  static  transmission  error  excitation, 
C=0.05,  Fm=0.1,  Qh=nT  and  a)  F=l.  b)  F=0.5. 
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higher  harmonics  of  Fj(t)  besides  the  fundamental  component  which  is  already 
included  in  the  analysis,  On  the  other  hand,  both  external  Fe(t)  and  internal  F,(t) 
excitations  may  exist  simultaneously.  These  two  cases  require  that  F(t)  in  equation 
(2,3)  be  rcfonnuluted  as  follows 

k 

F(‘)"Fm+  Fe(t)  +  Fj(t>=  Fm+  £  F  sin(jft  t+<>  ) 

j=  I  1 

+  £  (jQh)  FahjSUl(jQht +<|)hj)  (2.20) 

j  »  ] 

In  order  to  construct  analytical  solutions  of  this  equation,  we  must  investigate 
the  applicability  of  the  principle  of  superposition,  which  would  consider  first  each 
excitation  separately  and  then  superimpose  the  corresponding  responses  to  generate 
the  ovctull  frequency  response.  Two  cases  are  considered  and  the  analytical  results 
based  on  the  principle  of  superposition  ore  verified  through  digital  simulation. 

Consider  only  the  static  transmission  error  excitation  F4(t)  in  equation  (2.20) 

3  2 

with  3  harmonics,  i.e  Fc(t)  and  F(»)=Fm  +  £  (j^h)  Faj,jSin(jQht  +  4>hj)-  The 

H 

amplitudes  of  the  harmonics  are  selected  to  be  Fghi=0.05,  F4h2=0.02  and  F^sO.Ol, 
and  the  mean  force  l;m**0.1 .  First  each  excitation  component  is  considered  separately 
without  paying  any  attention  to  the  phase  angles  and  the  frequency  responses  are 
obtained  using  digital  simulation  us  shown  in  Figure  2.26.  Second,  all  there 
excitations  arc  included  simultaneously  and  the  ovcrull  response  is  compared  with  the 
previous  solutions  in  Figure  2,20,  The  ficquency  tesponsc  due  to  each  harmonic  of 


the  excitation  is  very  close  to  the  overall  solution  around  the  resonance  peaks 
governed  by  that  harmonic.  When  all  of  the  solutions  corresponding  to  each 
harmonic  are  added  algebraically  by  also  considering  the  phase  angles,  the  overall 
response  improves  as  shown  in  Figure  2.27.  This  response  is  in  the  form 
3 

qa(t)  =  L  tlajsm(^ht  +^hj)  "'here  qaj  is  the  alternating  displacement  when  only 

j=l 

the  j-th  harmonic  is  considered.  Figure  2.27  suggests  that  the  principle  of 
superposition  can  be  applied  to  a  gear  pair  with  backlash  provided  that  the  excitation 
frequencies  are  sufficiently  far  apart  from  each  other.  In  the  case  of  periodic  static 
transmission  excitation,  this  is  valid  as  all  excitation  frequencies  are  at  least  f2^  apart. 

Now,  the  principle  of  superposition,  which  has  been  already  verified  by  digital 
simulation,  can  be  used  to  obtam  the  approxunate  analytical  solutions  per  Section  2.4 
when  the  periodic  forcing  function  is  considered.  Figure  2.28  shows  the  analytical 
frequency  response  curve  for  Fm=0.1,  Fahj=0.05.  Fah2=0  02,  Fah?=0.01  and 
^=0.05.  Here  the  jump  discontinuity  is  seen  only  at  the  peak  Q()|  governed  by  Fahj 
since  Fm/Taj)?  and  Fm/Faj)^  arc  both  sufficiently  high  so  that  no  tooth  separation 
occurs  at  £2h2  and  respectively.  However,  the  jump  phenomenon  can  also  be 
seen  at  the  higher  harmonics  depending  on  the  force  ratios  Fm/Fahj,  j  *  1  and  £  in 
accordance  with  the  results  of  Section  2.6.  The  same  concept  can  be  applied  to  the 
superharmonic  components  of  the  external  torque  excitation  of  equation  (2.20) 
provided  Fallj=(). 

The  principle  of  superposition  is  now  extended  to  the  case  where  both  internal 
Fj(t)and  external  Fe(t)  sinusoidal  excitations  exist  simultaneously  as  given  by 
equation  (2.3)  provided  that  the  excitation  frequencies  i2()  and  f2-j-  are  not  close  to 
each  other.  In  a  real  geared  system,  i2(]  is  much  higher  than  iTj  which  implies  that 
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Figure  2. 20.  Comparison  of  the  fictjuency  rrsj>onse  due  to  fust  their  harmonics  of  the 

Matic  trancmission  error  for  the  cases  when  the  harmonic  *  ate  considered 
simultaneously  rutd  separately;  F(|(~(U,  f',1,1^0  0?,  J,,|1|,i^UU2  and 

Fjj^-O.Ol .  Ilicsc  results  arc  obtained  using  digital  simulation 
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Figure  2.27,  Comparison  of  digital  simulation  results  obtained  by  using  the  principle  of 
superposition  and  by  considering  all  three  harmonic  excitatiuns  applied 
simultaneously;  Fm<=0.1,  0.05,  ^,1*0.05,  Fth2li::0.02  and  F^sO.Ol.  In 

Uiis  case,  phase  angles  are  included  in  the  superposition. 
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Figure  2.28.  Comparison  of  the  frequency  response  qg  due  to  Fust  three  harmonics  of  the 
static  transmission  error.  This  curve  is  based  on  the  principle  of 
superposition  applied  to  the  analytical  solutions  which  have  been  constructed 
by  the  method  of  harmonic  balance,  Fm=0.1,  £=0.05,  F^jsO.05,  Fah2=0.02 

F3h3=0  01 
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Figure  2.29.  Frequency  response  qa  due  to  both  the  external  torque  and  the  static- 
transmission  error  excitations.  This  curve  is  based  on  the  principle  of 
superposition  applied  to  the  analytical  solutions  which  have  been  constructed 
by  the  method  of  harmonic  balance;  Fm=0.1,  ^=0.0.‘  ,  Fah=0.05.  Fal-=0.05, 
nh=2QT. 
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the  principle  of  superposition  should  be  suitable  for  this  case  Figure  2.29  illustrates 
the  frequency  response  solutions  due  to  the  fundamental  harmonic  component  of  both 
excitations  for  F|n=0. 1 ,  Fahj=FaTi=0.05,  q=0.05  and  by  assuming  that  Qh=2QT. 
Here,  the  jump  discontinuity  is  seen  at  two  frequencies.  However,  when  £2  h  =  £2  T 

or  when  one  of  the  superharmonic  peaks  of  Fj(t)  coincides  with  the  resonant  peak 
governed  by  Fe(t),  the  principle  of  superposition  will  no  longer  be  applicable.  In 
such  cases,  we  w  ill  use  the  digital  simulation  as  the  analytical  interpretation  of  these 
cases  is  yet  to  be  explored. 

2.9.  CONCLUSION 

This  analytical  study  on  the  non-linear  dynamics  of  a  spur  gear  pair  with 
backlash  as  excited  by  the  static  transmission  error  has  made  a  number  of 
contributions  to  the  state  of  the  art.  First,  difficulties  associated  with  the  digital 
simulation  technique  have  been  resolved  as  multiple  steady  state  solutions  at  a  given 
frequency  can  be  found  provided  the  entire  initial  condition  map  is  searched.  Second, 
new  frequency  response  solutions  for  the  gear  pair  have  been  constructed  using  the 
method  of  harmonic  balance.  Third,  our  mathematical  models  have  been  validated  as 
these  compare  well  with  two  previous  experimental  studies  and  the  key  parameters 
such  as  the  mean  load,  mean  io  alternating  force  ratio,  damping  and  backlash  have 
been  identified.  Forth,  the  chaotic  and  subharmonic  resonances  are  observed  if  the 
mean  load  is  too  small  for  a  lightly  damped  system.  Fifth,  mathematical  conditions 
for  tooth  separation  and  back  collision  have  been  established  which  are  compatible 
with  available  measured  data.  Sixth,  the  periodic  transmission  error  excitation  case 
has  been  analyzed  using  the  method  of  harmonic  balance  in  conjunction  with  the 
principle  of  superposition.  Finally,  on  a  more  fundamental  note,  our  study  enriches 


the  current  literature  on  the  clearance  non-linearity  or  vibro-impact  systems  as  the 
governing  non-linear  differential  equation  is  different  from  the  conventional  single 
degree  of  freedom  system  formulation. 
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CHAPTER  III 


NON-LINEAR  DYNAMIC  ANALYSIS  A  CEARED  ROTOR- 
BEARINC  SYSTEM  WITH  MULTIPLE  CLEARANCES 

3.1.  INTRODUCTION 

Mathematical  modeling  of  geared  rotor-bearing  systems,  being  an  essential  step  in 
designing  quiet  and  reliable  power  transmissions,  has  been  the  subject  of  numerous 
studies  over  the  past  few  decades.  Most  proposed  dynamic  models,  as  reviewed  by 
Ozguven  and  Houser  [14],  are  essentially  linear.  However  in  several  cases,  it  has  been 
experimentally  shown  that  the  geared  systems  exhibit  non-linear  behavior  [27-29,35]. 
For  instance,  vibro-impacts  are  observed  in  a  lightly  loaded  transmission  with  gear 
backlash  or  loose  bearings  [30,31,34],  Accordingly,  we  must  develop  non-linear 
mathematical  model  of  the  geared  system  -  this  is  main  focus  of  this  chapter  with 
emphasis  on  the  clearance  type  non-linearity  in  gears  and  rolling  element  bearings.  In 
Chapter  II  we  had  considered  the  single  degree  of  freedom  non-linear  model  of  a  spur 
gear  pair  with  backlash  and  investigated  the  effects  of  system  parameters  on  the 
vibrations  and  chaos  excited  by  the  static  transmission  error.  Applicability  of  the 
harmonic  balance  method  and  digital  simulation  technique  to  the  solution  of  the  steady 
state  response  has  been  demonstrated,  and  difficulties  associated  with  the  digital 
simulation  technique  when  applied  to  such  systems  governed  by  stiff  non-linear 
differential  equations  have  also  been  resolved  in  Chapter  II . 

Although  there  is  a  vast  body  of  literature  which  considers  a  single  degree  of 
freedom  system  with  clearances,  as  reviewed  earlier  in  Chapter  D,  studies  on  the  multi 
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degree  of  freedom  vibro-impact  systems  are  very  limited.  For  example,  Galhoud  et.al. 

(52)  considered  a  two  degree  of  freedom  translational  system  with  a  gap  and  found  the 
forced  harmonic  response  using  the  piecewtse  linear  technique.  Winter  and  Kojima 

[53]  also  used  the  same  technique  to  study  the  geared  systems  with  backlash 
However,  it  should  be  noted  ti.at  the  piecewise  linear  technique  can  not  predict  several 
non-linear  phenomena  such  as  subhaiir.onic  and  chaotic  responses  since  it  employs  the 
assumption  that  both  impact  and  no-impact  regimes  "are  repeated  in  an  identical  manner 
once  every  period  of  excitation"  [52]  Kucukay  [35]  developed  an  eight  degree  of 
freedom  model  of  a  helical  gear  pair  with  backlash  to  include  the  rocking  and  axial 
motions  of  rigid  shafts  and  radial  deflections  of  linear  bearings.  Lin  et.al.  [38]  included 
motor  and  load  inertias  in  a  three  degree  of  freedom  torsional  model.  Both  of  these 
studies  employed  the  digital  simulation  technique,  but  did  not  consider  a  number  of 
issues  of  primary  concern  in  non-linear  system  such  as  the  existence  of  multiple  steady 
state  solutions,  their  dependence  on  initial  conditions,  subharmonic,  quasi-periodic  and 
chaotic  responses,  etc.  A  few  of  these  ijsues  have  been  addressed  by  Singh  et.al  [34] 
on  the  gear  rattle  problem.  Comparin  and  Singh  [54]  have  also  used  the  digital 
simulation,  analog  simulation  and  harmonic  balance  method  to  analyze  coupled  impact 
pairs  assuming  that  the  modes  are  "weakly  coupled",  which  allows  the  system  to  be 
represented  by  a  combination  of  single  impact  pairs.  They  included  only  the  low 
frequency  external  torque  excitation  and  found  vne  steady  state  frequency  lesponse  at 
the  primary  resonance.  This  solution  has  then  been  used  to  analyze  the  neutral  gear 
problem  in  more  detail  [55 ] 

1  he  other  two  groups  of  studies  which  consider  multi-degree  of  freedom  systems 
with  continuous  non-linearities  [56-61]  and  periodic  excitations  due  to  mesh  stiffness 
variations  [29,62],  will  not  be  addressed  here  since  their  responses  are  significantly 
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different  from  the  clearance  non-lineai i;..s  as  discussed  by  Comparin  and  Singh  [33 J. 
The  literature  on  chaotic  vibrations  and  rolling  element  healings  will  be  discussed  later 

3.2.  PROBLEM  FORMULATION 
3.2.1.  Scope 

Figure  3.1  shows  the  generic  geared  system  considered  in  this  study.  It  consists 
of  a  spur  gear  pair  mounted  on  flexible  shafts  which  are  supposed  by  rollmg  element 
tarings  assembled  in  a  rigid  gear  box  Since  the  shafts  and  bearings  are  compliant, 
our  single  degree  of  freedom  model  of  spur  ^ear  pair  which  assumes  fixed  gear  centers 
is  obviously  not  suitable.  Instead,  a  three  decree  of  freedom  non-linear  model  as 
shown  in  Figure  3.2a  is  considered.  It  Includes  equivalent  stiffness  and  damping 
elements  representing  the  shaft  and  the  bearings.  The  corresponding  linear  model  has 
been  found  to  be  sufficiently  accurate  when  compared  with  a  fuhte  element  model  for 
eigen  solutions  piovided  that  the  gear  dynamic  response  (mesh  force,  dynamic 
transmission  error,  etc  )  is  of  major  concern  Gear  backlash  and  radial  clearances  in 
bearings  are  defined  analytically.  An  approximate  ncn-linear  bearing  model  is  also 
proposed.  Applicability  of  both  analytical  and  numerical  solution  technique*  to  ti  s 
problem  is  investigated  Several  key  issues  such  as  non-linear  modal  interactions  and 
differences  between  internal  static  transmission  error  excitation  and  external  torque 
excitation  are  discussed.  Parametric  studies  are  conducted  to  understand  the  effect  of 
system  pr  ameters  such  as  bearing  stiffness  kbj  to  mesh  stiffness  kh  ratio 

k  =  kbl  /  kh  ,  alternating  to  mean  force  ratio  Fa/Fm  and  radial  bearing  preload  to  mean 
force  ratio  Ft/Fni  on  the  non-linear  frequency  response.  A  criterion  is  used  to  classify 
ihe  steady  state  solutions,  and  the  conditions  for  chaotic,  quasi-periodic  and 
subhannonic  steady  state  solutions  are  detennined.  Tw’o  typical  routes  to  chaos 
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Generic  geared  rotor-bearing  system. 


observed  on  a  geared  system  are  also  identified.  Finally,  our  formulation  is  verified  by 
comparing  predictions  with  measurements  [27J. 

3.2.2.  Physical  Model  and  Assumptions 

The  three  degree  of  freedom  non-linear  model  of  the  geared  rotor  system  with 
gear  inertias  Ig|  and  Ig2,  gear  masses  mgl  and  mg2,  and  base  circle  diameters  dg|  and 
dg2,  as  shown  in  Figure  3.2a,  is  considered  here.  The  gear  mesh  is  described  by  a 
non-linear  displacement  function  fj,  and  viscous  damping  ch.  Friction  forces  at  the 
mesh  point  are  assumed  to  be  negligible  [18J.  Thus  the  transverse  vibrations  in  the 
pressure  line  direction  are  uncoupled  from  the  vibrations  in  the  direction  perpendicular 
to  the  pressure  line.  Bearings  and  the  shafts  that  support  the  gears  are  modeled  by 
equivalent  elements  with  viscous  damping  coefficients  c^  and  q,2  and  non-linear 
springs  defined  by  force-displacement  functions  fbl  and  fb2-  The  effects  of  the  prime 
mover  and  the  load  inertias  are  not  considered  assuming  that  these  inertial  elements  are 
connected  to  the  gear  box  through  soft  torsional  couplings.  Further,  we  assume  that 
the  system  is  symmetric  about  the  plane  of  the  gears  and  the  axial  motion  (parallel  to  the 
shafts)  is  negligible.  Like  the  spur  gear  pair  model  of  Chapter  II,  both  low  frequency 
external  excitation  due  to  torque  fluctuations  and  high  frequency  internal  excitations  due 
to  the  static  transmission  error  e(t)  are  considered  in  the  formulation.  Input  torque 

fluctuation  is  included,  but  the  output  torque  is  assumed  to  be  constant,  i.e. 
Tgi(t)  =  Tgim  +  Tgla(t)  and  Tg2(t)  =  Tg2m.  External  radial  preloads  F^  and  Fj,2  are 

also  applied  to  both  rolling  element  bearings. 
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3,2  Non-linear  models  of  the  geared  rot  or -bearing  system,  a)  tluee  degree  of 
freedom  model,  b)  two  degree  of  freedom  model 
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p(t)  =  x(t)  +  ygl(t)-yg2(t)-e(t) ; 


(3.2a) 
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A  dimensionless  form  of  equation  (3.2b)  is  obtained  by  letting 
ygl(t)  =  ygi(t)/bc,  p(0  =  p(t)/bc,  o>n  =  yj kb  /mcl  ,  cobi  =  ^kbi  Jm~  (i=l,2)  and 

t  =  o)nt  where  bc  is  the  characteristic  length.  Here,  we  consider  harmonic  excitation  for 
both  e(t)  and  FaT(t)  as  e(t)  =  esin(Qbt  +  <|>h) ,  FaT(t)=:  FaTsin(flyt  +  <j>T)  where  Qb 

and  Qj  are  the  fundamental  excitation  frequencies  of  internal  displacement  and  external 
torque  fluctuations,  respectively.  Further,  define  dimensionless  excitation  frequencies 
f/h  ~  /  (On  and  =  £2x  /®n  to  yield  the  following  dimensionless  governing 

equation  of  motion. 
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{F(t)}  =  {F}m  +  {F(t)}4  +  {F(i)}c 
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where  (i*l  ,2)  and  Fm,  are  the  dimensionless  components  of  the  mean  force  vector 
(F)m,  and  Fay  and  pertain  to  the  alternating  external  excitation  (F(t))r  and  internal 
excitation  |F(t)}j  force  vectors,  reipcctivcly. 

3.2.4.  Modeling  of  Non-lineurilies 

The  non-linear  displacement  functions  ft,|(ygi).  !■!,  2  and  fj,(p)  in  equation 
(3.3a),  which  represent  the  bearing  radial  and  gear  mesh  stiffnesses,  respectively, 
should  be  defined  explicitly  before  solving  the  non-linear  equations.  Here,  f|,(p)  is, 
defined  as  a  clearance  type  dead  space  function  with  a  backlash  2bj,  and  linear  time- 
invariant  mesh  stiffness  of  one,  in  the  dimensionless  form 
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(3.4) 
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For  the  i-th  rolling  element  bearing,  radial  force  Pyi  versus  displacement  y g 
relationship  under  the  static  loading  condition  is  defined  as  [63-65] 
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(3.5) 


where  k,,  is  the  inner  contact  stiffness,  ar  is  the  angular  position  of  the  r-th  rolling 
element  in  contact,  2bj,|  is  the  radial  clearance  of  the  i-th  bearing,  n  is  the  power  of  the 
non-linear  force  displacement  relationship  (n=1.5  for  ball  bearings  and  n=10/9  for 
roller  bearings)  and  11  is  the  total  number  of  rolling  elements  in  contact  under  loaded 
conditions.  Now  the  dimensionless  bearing  displacement  function  fbi(ygj)  of  equation 
(3,3a)  is  obtained  from  equation  (3.5)  as 
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fbi(y*i)  =  7-^r  =  i 
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Figure  3.3  shows  the  non-iinear  function  fgi(ygi)  for  a  roller  bearing  with  n=10/9, 
kt=lxl08  N/mlO/9,  total  number  of  rollers  Z=15  and  ^=0.01  mm.  In  Figure  3.3,  we 
note  almost  a  linear  relationship  for  large  displacements,  say  ygi>3bbi/bc.  Since  the 
degree  of  non-linearity  is  not  very  significant,  equation  (3.6)  can  be  approximated  by  a 
piece-wise  linear  function,  similar  to  fh(p)  given  by  equation  (3.4),  in  order  to  simplify 
the  analysis  considerably.  Figure  3.3  also  illustrates  two  linear  approximations  A  and 
B  beyond  the  clearance  for  fbi(ygi)  in  the  form 
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Note  that  both  approximations  A  and  B  differ  in  clearances  bj,iA  and  b^jg  but  have  the 
same  slope  as  the  exact  bearing  stiffness  curve  for  ygj>3btjj/bc.  i=l  ,2,  which  is  unity  in 
the  dimensionless  form.  The  validity  of  these  approximations  will  be  given  later. 
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Figure  3.3 .  Exact  and  approximate  bearing  deflection  functions. 
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2.3.  CORRESPONDING  LINEAR  MODEL 

As  a  limiting  case,  equations  of  motion  of  the  corresponding  linear  system  are 
obtained  by  substituting  fbi (ygi )  =  ygi  .  for  i=l  and  2,  and  fb(p)  =  p  in  equation  (3.3a) 

'i  o  o][ygi(Oj  r;n  o  Ci3lfygi(o 

o  1  0  jyg2(0*+2  0  C22  —  C23  ■  yg2(0  ■ 

--1  1  1  P(0  .0  0  C33.  p(0 

'* it  o  ^13'fygi(t) 

+  0  K22  -k23  •  yg2(t)  ={F(t)}  (3.8a) 

0  K33  J  p(t) 

or  in  the  matrix  form  with  { q(t) }  as  the  displacement  vector 

[MJ(q(t)|  +[C]|q(t)J  +  [K]|q(t)|  =  |F(t)}  (3.8b) 

where  the  mass  [M],  damping  (C]  and  stiffness  [K]  matrices  are  all  positive  definite. 
These  matrices  are  asymmetric  due  to  ygj  and  yg2  terms  in  the  last  row  of  matri? 

equation  '3.8a).  For  this  linear  system,  the  mean  and  alternating  components  of  the 
motion  can  be  separated  by  letting  ygi(t)=ygjm+ygia(t),  i=l,2  and  p(t)=pm+pa(t). 
Hence,  equation  (3.8a)  is  rearranged  in  terms  of  the  alternating  motion  as 

[M]{qa(t»  +  [C](qa(t))  +[K]|qa(t)|  =  (F(t))j  +  (F(t))e  (3.9) 

The  natural  frequencies  cor  and  the  modes  {yr}  are  calculated  by  considering  the 
corresponding  eigen-value  problem.  The  forced  harmonic  vibration  response  is  then 
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obtained  by  the  modal  expansion  technique  in  the  following  form 


<qa(t» 


3 

I 


r=l 


(¥fH¥r>T 

(a)J-^)  +  jahcr 


Fah^h  sin(Qht  +  4>h) 


+  Z  — __  I  ¥r>l  ¥  r ) - •  FgjQy  s*r>(^Tt  j  =  V-T  (3.10a) 

r=t  (tor  -QT)  +  jQTcr 


cr=^-l¥r)T[C)l¥5);  5re  = 


'rs 


1;  r  =  s 
0;  r^s 


(3.10b) 


Table  3.1  shows  the  natural  frequencies  of  the  three  degree  of  freedom  linear 
system  given  by  equation  (3.9)  and  the  ones  found  by  finite  element  method  for  three 
different  k  =  kbi/kh  values  with  mgj=l  kg,  Igj=0.0008  kg/m2,  dgi=0.08  m. 

Cjj=0.01,  ^3=0.0125,  ^33=0.05,  i=l,  2  and  kh=2xl08  N/m.  As  evident  from  Table 
3.1,  the  three  degree  of  freedom  and  finite  element  models  [15]  result  in  virtually  the 
same  natural  frequencies.  This  concludes  that  the  three  degree  of  freedom  model  is 
indeed  suitable  for  the  geared-rotor  bearing  system.  The  first  \j/j  and  third  \j/in  natural 
modes  are  coupled  transverse -torsional  modes  while  ¥y  *s  purely  transverse  type  [15]. 
The  second  natural  mode  \|/jj  is  not  excited  by  e(t)  in  this  particular  case,  since  the  gear 
ratio  .Vg=dg2/dgi  is  one  for  a  symmetric  pair.  Therefore  only  two  peaks  should  exist  in 
the  frequency  response  spectrum  which  will  be  presented  in  Section  3.5. 
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Table  3.1.  Natural  frequencies  of  the  corresponding  three  degree  of  freedom  (3- 


DOF)  linear  model . 


k=kbi/kh 


1 

5 

10 

Vr 

u>r 

3-DOF/FEMt 

3-DOF/FEM 

3-DOF/FEM 

first  transverse- 
torsional  coupled 

o>l 

0.396/0.402 

0.500/0.512 

1.262/1.275 

purely  transverse 

con 

0.758/0.795 

1.118/1.157 

1.475/1.560 

second  transverse- 
torsional  coupled 

<°in 

0.880/0.930 

1.581/1.692 

1.796/1.868 

t  Finite  Element  Method  (FEM)  from  Reference  [15]. 
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3.4.  TWO  DEGREE  OF  FREEDOM  SYSTEM  STUDIES 
3.4.1.  Equations  of  Motion 

As  a  first  example,  we  reduce  the  three  degree  of  freedom  transverse-torsional 
semi-definite  model  to  a  two  degree  of  freedom  non-linear  model.  This  is  obtained  by 
clamping  one  of  the  gear  centers  as  shown  in  Figure  3.2b,  i.e.  one  of  the  transverse 
displacements  is  assumed  to  be  7.ero,  say  ygj=0  which  physically  corresponds  to  a 
system  with  one  gear  gl  mounted  on  a  rigid  shaft  which  is  supported  by  very  stiff 
bearings  while  the  other  gear  g2  is  assembled  with  compliant  lumped  shaft  and 
bearings.  Equations  of  motion  of  the  reduced  order  system  are  obtained  from  equation 
(3.3a)  by  letting  ygl=0. 


1 

1 


yg2«)j 

p(t)  J 


~  C23 

C33  p(0  J 


|fb2(yg2)| 
1  fh(p)  J 


(3.11) 


The  intent  here  is  to  simplify  the  physical  system  in  ordei  to  investigate  several 
key  issues  in  depth.  Specifically,  the  objectives  arc  to:  a)  justify  using  approximate 
bearing  models  proposed  in  Section  3.2.4  by  comparing  the  frequency  responses 
excited  by  e(t ) ,  b)  show  applicability  of  the  harmonic  balance  method  (KBM)  to  solve 
non-linear  system  equations  and  compare  its  predictions  with  the  results  yielded  by  the 
digital  simulation  technique,  c)  study  interactions  between  system  non-linearities  for 
both  weakly  and  strongly  coupled  modes,  and  d)  compare  the  steady  state  frequency 
response  spectra  due  to  internal  static  transmission  error  excitation  |F(t))j  and  external 
torque  excitation  (F(t))e, 


1 10 


3.4.2.  Solution 

Both  analytical  and  numerical  solution  techniques  which  we  h*ve  used 
successfully  for  the  spur  gear  pair  problem  are  employed  again  to  solve  equation 
(3.11).  First,  approximate  analytical  solutions  are  constructed  using  the  harmonic 
balance  method  [33,51]  with  only  the  bearing  non-linearity  asstuning  no  gear  backlash, 
i.e.  1^=0  or  fh(p)=P  For  the  harmo.uc  excitation  given  by  equation  (3.3) ,  we  assume 

p(t )  =  Pm  Pa  sin(f2|)t  +  <|>p )  i  yg2(0  =  yg2m  +yg2asill(ftht  +  4>g2)  (3.12a, b) 

where  subscripts  m  and  a  represent  mean  and  alternating  components  of  the  steady  state 
response,  and  <bp  and  <t*g2  are  the  phase  angles.  The  non-linear  bearing  function  is 
approximated  as 

^b2^g2 )  =  Nmyg2m  +  ^ayg2a  sin(Q^t  +  $g2  )■  (3.13) 

where  Nm=Nm(yg2m,yg2a)  and  Na=Na(yg2m-yg2a)  are  the  describing  functions  which 
are  given  in  Chapter  II;  these  need  to  be  defined  for  ea>  h  impact  regime  (no  impact, 
single  sided  impact,  double  sided  impact).  Governing  frequency  response  equations 
are  obtained  by  substituting  equations  (3.12)  and  (3.13)  into  equation  (3.11)  and 
equating  the  coefficients  of  like  harmonics: 

[A]-f(2;22Qhrl«hFah 

Pa  = - He - ;  t3.14a> 

V  A 


111 


A-l-QjtajAj  +4^hC2:C;3  +  a2A1  ^  4fthC22A2)<L 


+(-2P.^22K23  -2^h^23Ai  +  2f2h^3^  +8^1^22^33^  i  (3.14b) 


A1  =K22Na  A2  =  1  ~^h  i  Pm  =  Fm  ; 


(3. 14c,d,e) 


yg2a  = 


*23  ':'(2^bC23)2 


L( k23A1  “ 4ftb^22^23 )  +  tf°hC22K23  +  2tib^23Al )  J 


1/2 


(3.140 


yg2m  - 


^b2  *23Pm 

*22Nm 


(3.1 4g) 


The  overall  frequency  response  is  obtained  by  solving  equation  (3.14)  numerically  for 
each  impact  regime  separately. 

Second,  the  digital  simulation  technique  is  used  to  solve  equation  (3.1 1)  for  no 
gear  backlash.  A  5th-6th  order  variable  step  size  Runge  Kutta  numerical  integration 
algorithm  (50],  which  has  already  been  employed  successfully  in  solving  similar 
problems,  is  used  here.  Figure  3.4  compares  frequency  response  spectra  obtained 
using  the  harmonic  balance  method  and  digital  simulation.  Although  both  methods  are 
in  a  very  close  agreement  for  the  case  considered  in  Figure  3.4,  one  should  be  aware  of 
the  following:  (i)  several  problems  may  exist  in  the  application  of  the  numerical 
integration  method  to  the  clearance  type  problems,  hence  caution  must  be  exercised, 
and  (ii)  the  harmonic  balance  method  is  incapable  of  predicting  chaotic  and 
subharmonic  responses.  Accordingly  for  further  analysis,  we  will  use  only  the  digital 
simulation  technique. 
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amplitude  amplitude 


a) 
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'"c"'  digital  simulation 
A  HBM 


frequency 


Figure  3.4.  Comparison  of  harmonic  balance  method  and  digital  simulation  technique. 

Fm  =  l-0,  Fah=0.5,  Fb2=0,  bb=0,  bb2  =  bc,  ;33=0.05.  ;23=0.0125, 
^22=0  01,  k  =  1;  a)  yg2a  versus  fth,  b)pa  versus  flh. 


3.4.3,  ViiliUultun  of  Non *1  Incu r  ItcurlnK  Model 

We  examine  (he  validity  of  two  approximations  for  fb|(yg|).  «*>  auggeited  in 
Section  3,2  4  and  iiliown  in  Figure  3,3,  by  assuming  hh*K)  am)  F^wO,  Equation 
(3.11)  ia  solved  using  the  digital  simulation  technique,  given  l'mt.|  ,0,  F^O.O.V 

Fb2*°-  ^22"°  ,01, 4))"0,05,  £21^01 23  and  K22s,*2j*,0,23,  for  three  hearing  model! 
(exact,  approximations  A  approximation  B;  »l  ■63h1) 

Figures  3.3a  and  3  5b  compare  versus  il(l  and  p,  versus  U|,  e|*tctra  for  aJI  thiee 
hearing  models.  Both  approximations  me  in  a  very  close  agrer/meni  with  the  exact 
hewing  model  over  the  entire  frequency  range  ihe  frequency  and  the  amplitude  of  liie 
jump  at  the  first  peak  is  predicted  accurately  by  hulls  approximate  models  Hence  we 
conclude  tiiai  the  approximate  piecewise  linear  model  can  I*  used  for  a  Iteanng  without 
losing  any  accuracy  and  the  dynamic  behavior  is  not  sensitive  to  the  minor  variations  in 
hearing  clearances  ty,, 

3.4.4,  Nondliieiir  Modul  litter  uctlons 

Although  both  natural  modes  of  the  corresponding  linear  model  are  transverse- 
torsional  coupled  modes,  first  mode  q/j  Is  more  etc  pendent  on  transveise  vibrations  than 

the  second  peak  Accordingly,  it  is  expected  that  (tearing  nun  linrarliiea  should  affect 
the  first  peak  directly  while  gear  backlash  should  dominate  the  second  peak  in  Figure 

3.4,  Hie  validity  of  this  claim  depends  on  tire  nature  of  non  linearities  and  system 
parameter  values,  as  illustrated  below. 

First,  a  set  of  parameters  is  selected  so  that  natural  frequencies  «>j  and  orj|  of  tire 
corresponding  linear  system  are  far  away  from  each  other;  osj *^0  4^7  anti  ii>||«l  144 
for  K22*0,23,  h23«*0,23,  4;;*»0,0I.  4;p*(UJl23,  and  03.  Frecpieiicy  res|»oii!ie 
characteristics  of  die  sysicm  widi  a  forcing  function,  F(l  ■*!  .0  and  Fg|,*(),3  are  evaluated 


i  M 


frequency 


Comparison  of  frequency  responses  obtained  by  using  the  exact  and 
approximate  bearing  displacement  functions  given  in  Figure  3.3; 

Fah=0.5.  Fb2=0-  hh=0,  bb2=bc,  ;33=0.05,  ;23=O  0125, 
C22r()  0 1  .  k  -•  1 ;  a)  yg2a  versus  £2h,  b)  pa  versus  Qb. 


for  the  following  three  cases:  a)  no  gear  backlash,  1^=0,  non-linear  bearing,  b^2=bc;  b) 
gear  backlash,  bj,=  10bc,  linear  bearing,  1^2=0;  and  c)  gear  backlash  and  non-linear 
bearing  with  b|1=10b|)2=10bc.  In  the  last  two  cases,  gear  backlash  is  taken  to  be  10 
times  the  bearing  clearance  as  it  is  a  typical  order  of  magnitude  value  for  many  geared 
rotor  systems  supported  by  roller  bearings.  Figures  3.6  and  3.7  show  the  frequency 
responses,  yg2a  versus  and  pa  versus  respectively,  for  all  three  cases  given 
above;  also  superimposed  is  the  linear  frequency  response  curve  for  comparison 
purposes.  As  shown  in  Figure  3.6a  and  3.7a,  a  jump  discontinuity  is  observed  at  the 
first  peak  only  for  case  a.  In  the  case  with  only  gear  backlash  non-linearity,  jump  is 
seen  at  the  second  peak  while  the  first  peak  is  continuous  as  shown  in  Figures  3.6b  and 
3.7b;  this  supports  the  claim  of  torsional  mode  dominance  of  the  second  peak.  Finally, 
when  both  non-linearities  are  included  simultaneously,  both  peaks  become  non-linear 
as  shown  in  Figures  3.6c  and  3.7c.  For  such  a  system,  natural  modes  are  essentially 
"uncoupled"  and  the  interactions  between  component  non-linearities  are  negligible.  In 
summary,  for  case  a  of  Figures  3.6a  and  3.7a,  a  linear  model  can  b  u*en  ;o  predict  the 
frequency  response  beyond  (Oj  ,  say  Q^>  0.6.  Similarly,  a  linear  model  could  be 
suitable  for  case  b  for  Qh<0.8  as  shown  in  Figuies  3.6b  and  3.7c.  But  for  case  c,  a 
linear  model  is  suitable  only  far  away  from  the  resonances. 

Second,  choose  a  dataset,  say  K22=0.64,  <23=0.10,  ^2=0. 01,  ^23=0.01  25,  and 
^33=0.05,  such  that  two  natural  frequencies  for  the  corresponding  linear  system  are 
brought  closer  to  each  other:  o>j=0.720,  o>j|=l  .101 .  We  again  consider  the  same  cases 
a,  b  and  c  and  compare  the  frequency  response  for  Fm-2.0  and  F^sO.5  in  Figures  3.8 
and  3.9.  Here,  each  non-linearity  affects  both  modes,  resulting  in  jumps  at  each  peak. 
Therefore,  for  this  set  of  parameters,  modes  are  considered  to  be  "coupled"  because  of 
dynamic  interactions.  Accordingly  both  non-linearities  must  be  included 
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b) 
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frequency 


Figure  3.6.  yg2«  v*i’*us  U|,  plots  foi  the  case  of  uncoupled  modes;  a)  no  backlash, 
non-linear  bearings,  b)  backlash,  linear  bearings,  c)  backlash,  non-linear 
bearings; - corresponding  linear  response 
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5.0 


b) 


Figure  3,7.  pa  versus  f2h  plots  for  the  case  of  uncoupled  modes,  a)  no  backlash,  non¬ 
linear  bearings,  b)  backlash,  linear  bearings,  c)  backlash,  non-linear 
bearings; - corresponding  linear  response. 
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frequency 


Figure  3.8.  yg2a  versus  Qh  plots  for  the  case  of  coupled  modes;  a)  no  backlash,  non¬ 
linear  bearings,  b)  backlash,  linear  bearings,  c)  backlash,  non-linear 
bearings;  corresponding  linear  response. 
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frequency 


Figure  3.9.  pa  versus  £2h  plots  for  the  case  of  coupled  modes;  a)  no  backlash,  non¬ 
linear  bearings,  b)  backlash,  linear  bearings,  c)  backlash,  non-linear 
bearings; - corresponding  linear  response. 


simultaneously  in  the  dynamic  model. 


3.4.5.  Internal  versus  External  Excitation 

Here  we  will  apply  only  one  type  of  excitation  at  a  time  and  compare  frequency 
response  characteristics  of  the  system  due  to  sinusoidal  excitations  {F(t))j  and  { F(t))e 
in  equation  (3.3a).  Here,  only  the  bearing  non-linearity  will  be  considered  assuming 
no  tooth  separation,  b},=0.  In  the  case  of  internal  excitation,  amplitude  of  alternating 
component  is  dependent  on  ftb .  Conversely,  the  amplitude  of  alternating  external  force 
is  frequency  independent.  Therefore,  for  a  constant  force  ratio  Fai/Frn,  i=h,T,  the 
overall  amplitude  ratio  ^hFjjf/Fn,  for  internal  excitation  varies  with  excitation  frequency 
wheieas  the  overall  amplitude  ratio  FaT/Fm  for  external  excitation  remains  constant 
similar  to  the  single  degree  of  freedom  model  of  Chapter  II. 

Figures  3.10a  and  3.10b  compare  yg2a  versus  ft  and  pa  versus  ft  spectra  for 
ftt=fth=ft,  Fm=1.0,  FaT=Fah=0.5,  Fb2=0,  C22=0.01,  C23=0.0125,  C33=0.05  and 
k22=k23=0-25.  As  shown  Figures  3.10a  and  3.10b,  both  excitations  yield  the  same 
values  at  ft=1.0  since  Fm/FaT=Fm/ft2Fab  at  ft=1.0.  For  ft>l  the  internal  excitation 
gives  larger  amplitudes  than  the  external  excitation,  but  the  converse  is  true  for  ftcl . 
Another  important  difference  is  that  for  increasing  ft,  there  are  two  jumps  (  a  jump-up 
from  the  no-impact  regime  to  the  double  sided  impact  regime  at  ft=0.3  and  a  jump- 
down  from  the  double  sided  to  the  single  sided  impact  regimes  at  ft=0.4)  for  the 
external  excitation;  but  the  double  sided  impact  regime  is  not  seen  in  the  internal 
excitation  case  and  only  one  jump  from  the  no-impact  to  the  single  sided  impact  is 
found  at  ft=0.35. 
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amplitude  amplitude 


a) 


Figure  3.10.  Comparison  of  frequency  responses  due  to  internal  static  transmission 
error  and  external  torque  excitation;  Fm=1.0,  Fah=0.5,  Fb2=0,  bh=0, 

bb2=bc*  C33=0  05>  C23=00125,  ^22=°  k  =  1;  a)  y^2a  versus  fth,  b) 
pa  versus  fib. 
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3.5.  THREE  DECREE  OF  FREEDOM  SYSTEM  STUDIES 
3.5.1.  Classification  of  Steady  State  Solutions 

The  steady  state  response  results  excited  by  internal  force  { F(t) )  j  at  each 
frequency  Qj,  have  been  categorized  into  the  following  five  groups:  (i)  harmonic  01 
nearly  harmonic  solution  at  the  same  period  tp=2rt/Oh  as  that  of  forcing  function 
(period-one),  (ii)  non-harmonic  period-one  solution,  (iii)  subharmonic  solution  with 
period  ntp,  n>l,  (iv)  quasi-periodic  solution,  and  (v)  chaotic  solution  (non-periodic, 
n  — ►  °o).  The  solution  classification  criteria  are  based  on  time  histories,  phase  plane 
plots,  Poincare  maps  and  Fourier  spectra  {66-71].  Figures  3.11  through  3.15  illustrate 
different  types  of  steady  state  solutions  obtained  from  the  three  degree  of  freedom  non¬ 
linear  model  of  Figure  3.2a.  As  shown  in  Figure  3.11,  period-one,  non-harmonic 
solutions  have  a  non-circular  (non-elliptic)  phase  plane  plot  and  repeat  themselves  at  tp. 
Figure  3.12  shows  the  time  histories,  phase  plane  plots  and  Poincare  maps  for  a  2tp 
subharmonic  solution.  In  this  case,  the  Poincare  map  consists  of  two  discrete  points. 
Similarly,  the  6tp  (ultra)subharmonic  response  has  six  points  on  the  Poincare  map  as 
shown  in  Figure  3.13.  Figure  3.14  illustrates  a  typical  quasi-periodic  response 
(combination  oscillations)  which  consists  of  two  or  more  "incommensurate" 
frequencies  [66,67],  Quasi-periodic  solutions  result  in  closed  orbits  on  the  Poincare 
map  as  shown  in  Figure  3.14c.  The  chaotic  responses  which  are  defined  by  a  non¬ 
periodic  time  history  and  as  many  points  as  the  number  of  cycles  considered  in  the 
analysis  on  the  Poincare  map  are  illustrated  in  Figure  3.15.  Figures  3.16a-e  show  the 
FFT  spectra  corresponding  to  the  time  histories  given  in  Figures  3.1 1  to  3.15.  For 
period-one  non-harmonic  solution  of  Figure  3.1 1,  corresponding  spectrum  has  peaks  at 
mo)  where  to  is  the  fundamental  frequency  and  m  is  a  positive  integer,  as  shown  in 
Figure  3.16a.  The  ntp  subhamionic  solutions  contain  peaks  at  the  frequencies  mw/n. 
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For  instance,  spectra  of  the  2tp  and  6tp  solutions  of  Figures  3.12  and  3.13  have  peaks 
at  frequencies  mto/2  and  mto/6  respectively  in  Figures  3.16b  and  3.17a.  The  quasi- 
periodic  solution  of  Figure  3.14  consists  of  two  fundamental  frequencies  (Oj  and  Wj  at 
a  ratio  /  (&2  =8  (an  irrational  number)  and  there  are  peaks  at  the  combination 
frequencies  nuoi+ro^,  m,r  =  0,  ±  1,±  2....  Finally,  a  characteristic  broad  band 
spectrum  is  obtained  in  Figure  3.18  when  the  solution  as  shown  in  Figure  3.15  is 
chaotic. 

3.5.2.  Routes  to  Chaos 

In  Chapter  D,  we  had  investigated  the  effect  of  mean  load  Fm  and  damping  £  on 
chaotic  response  of  a  gear  pair.  It  has  been  shown  that  chaos  typically  exists  in 
lightly  damped  and  lightly  loaded  gear  pairs.  Now,  two  different  routes  to  chaos  for 
the  three  degree  of  freedom  system  are  illustrated. 

A.  Period-Doubling  Route  to  Chaos :  This  consists  of  a  sequence  of  bifurcations 
of  the  periodic  response  to  another  periodic  response  with  twice  the  period  of  the 
original  response  due  to  a  change  in  one  of  the  system  parameters  (66,  68-70].  Figures 
3.19  to  3.23  demonstrates  this  with  a  change  in  excitation  frequency  Cth  given  Fm=0.1, 
Fah=0-05,  Fj,,=0  i=l,2,  Cj  1  =  ^22=0.01 ,  ^3=^23=0.0125,  £33=0.05,  Kj j=K22=1 -25, 
Ki3=K23=0.25,  bc=bjj  and  bt>j=0,  i=l,2.  At  0^=1.500,  with  zero  initial  conditions, 
we  note  a  tp  harmonic  solution  in  Figure  3.19.  This  tp  solution  bifurcates  to  2tp 
subharmonic  solution  in  Figure  3.20  when  f2b  is  reduced  to  £2^=1 .48.  Furthennore  at 
£2h=1.44  a  4tp  solution  is  obtained  which  is  further  transformed  to  a  8tp  response  at 
£lh=1.402  as  shown  in  Figures  3.21  and  3.22  respectively.  At  Qh=l  4,  however  the 
steady  state  solution  becomes  chaotic  as  evident  from  Figure  3.23.  Also  note  from 
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bb2=bc,  C33=0.°5,  ^13=^23=0  0125.  Cl  !=C22=0  01.  k=5,  Oh=0.6;  a) 


time  histories,  b)  phase  plane  plots. 


Figure  3.12.  A  period  - 
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Figure  3.17,  FFT  spectra  of  y^U)  for  the  responses  given  in  Figures  3,13  to  3,14  a) 
period-six,  b)  quasi-llnear, 
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frequency 


Figure  3.18.  FFT  spectra  ofyg2(t)  for  the  chaotic  response  given  in  Figure  3.15 
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spectra  of  Figures  3.19  to  3.23  that  the  single  peak  spectrum  is  transformed  to  a  broad 
band  spectrum  corresponding  to  chaos  through  period-doubling. 

B.  Quasi-periodic  Route  to  Chaos:  In  this  case,  instead  of  going  through  the 
period  doubling,  the  response  undergoes  a  sequence  of  Hopf  bifurcations  with  a 
change  in  a  system  parameter.  First  a  quasi-periodic  response  and  finally  a  chaotic 
response  are  obtained  [66,71];  Figures  3.24-3.27  shows  such  a  route  to  chaos.  3tp 
solution  at  1 . 1  bifurcates  to  a  quasi-periodic  solution  when  is  increased  to  1.2 
as  shown  in  Figure  3.24  and  3.25.  A  further  increase  in  to  1 .3  yields  a  deformed 
closed  orbit  Poincare  map  as  shown  in  Figure  3.26.  It  is  then  transformed  to  a  chaotic 
strange  attractor  at  £lh=l .4  in  Figure  3.27.  Similar  observations  are  evident  from  the 
Fourier  spectra  shown  in  Figures  3.25-3.27;  the  spectrum  consists  of  combination 
frequencies  for  the  quasi-periodic  solution,  but  it  finally  changes  to  a  broad  band 
spectrum  characteristic  of  chaos. 

3.5.3.  Parametric  Studies 

The  three  degree  of  freedom  non-linear  model  of  Figure  3.2a  is  used  to  study  the 
effect  of  several  system  parameters  such  as  the  stiffness  ratio,  k  =  kbi  /kh  i  =  1,2, 

radial  bearing  preload  to  mean  load  ratio,  Fbj/Fm  (i=l  ,2)  and  alternating  load  to  mean 
load  ratio  Fjb/F,,,  on  the  non-linear  dynamic  behavior.  A  geared  rotor-bearing  system 
of  one-to-one  gear  ratio  (vg=l)  with  mgj=l  kg,  Igj=0.0008  kg/m2,  dg,=0.08  m,  i=l  ,2, 
and  mesh  stiffness  kjj=2xl08  N/m  is  selected.  Five  percent  mesh  damping  and  one 
percent  bearing  damping  values,  which  are  somewhat  realistic  [72,73],  are  used.  i.e. 
Cj  1=^22=0.01,  ^13=^23=0.0125  and  £33=0. 05.  Both  non-linearities  are  considered 
separately  in  this  parametric  study.  The  case  when  both  bearing  non-linearities  and 
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Figure  3.22.  Frequency  doubling  route  to  chaos;  Fm=0. 1,  Fafo=0.05,  F^O,  bh=bc, 
bb2~0 ,C33=005,  ^13=^23=0  0125,  Cll=^22=0  01.  k=5;  period-eight 
response  at  ty,=  1.402. 


Figure  3.24.  Quasi-periodic  route  to  chaos;  Fm=0. 1.  Fah-0.05.  Ff,2=0.  hh=bc=0.1. 

bb2=0  ^33=0.05,  ^i3=C23=0.0125,  Cll-C22=^  ^1,  ^=1;  Pen°d  three 
response  at  fih=  1.1. 


bb2=0,  ^33=0.05,  ^i3=C:?=0  0125.  Cn=C22=0  01.  k  =  l ;  quasi-pcriodic 
response  at  ^=1 .2. 
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Figure  3.26.  Quasi-periodic  route  to  chaos;  Fm=0.1,  Faj,=0.05.  F'b2=0.  ht)=bc=0.l. 

bb2=0,  ^33«0.05,  Ci?=^2*=(,°125-  9ll=^22=^  01*  ^=1;  iransition  from 
quasi-periodic  to  chaotic  response  at  fth=i.3. 
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Figure  3.27.  Quasi-periodic  route  to  chaos;  Fm=0  1.  Fah=0.05.  Ft,2=0,  b|,=bc=0.1, 
bb2=0(  C?3=0-05,  Cn=C.23=0.01 25.  |  =  ^2 2  =0  0 1 .  k=l;  chaos  at 


Qh=1.4. 
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backlash  exist  simultaneously  is  not  considered  here  as  it  has  already  been  examined  in 
Section  3.4  using  the  two  degree  of  freedom  model. 

A.  Gear  Backlash.  Linear  Bearings:  For  the  non-linear  model  of  Figure  3.2a  with 
linear  bearings  (bjjj=0),  bj,  is  used  as  the  characteristic  length  bc.  Figures  3.28a  and 
3.29a  show  ygla(£2b)and  pa(Qb)  spectra  respectively  for  the  soft  bearings  (k  =  1)  with 
Fm=0.1  and  Fab/Fm=0.5.  Here  yg2a  spectrum  is  not  included  since  ygia=yg2a  (ygi(0=* 
yg2(t))  for  a  one-to-one  gear  ratio.  As  evident  from  Figures  3.28a  and  3.29a,  yg]a  at 
the  Fust  peak  is  larger  than  that  seen  at  the  second  peak.  Conversely,  pa  is  much  lower 
at  the  first  peak,  implying  that  first  mode  is  dominated  by  the  transverse  vibrations 

whereas  the  second  mode  is  dictated  by  the  torsional  vibrations  for  soft  bearings.  This 
corresponds  to  a  large  dynamic  bearing  force  Fdbi  =*  kbiygj .  but  a  small  dynamic 

transmission  error  (p(t)+e(t))  around  the  fust  natural  frequency.  Jump  phenomenon  is 
observed  only  at  the  second  peak  which  is  governed  by  the  gear  mesh  where  the 

A 

excitation  e(t)  is  applied.  This  indicates  that  for  k  =  1  the  modes  are  weakly  coupled. 
Accordingly  the  gear  mesh  non-linearity,  which  forces  the  second  peak  to  be  non¬ 
linear,  does  not  affect  the  linear  characteristics  of  the  Fust  peak.  However  the  modal 
coupling  becomes  stronger  with  an  increase  in  k  as  evident  from  the  jumps  seen  at  both 
peaks  in  Figures  3.28b  and  3.29b  for  k  =  5  and  Figures  3.28c  and  3.29c  for  k  =  10. 
Torsional  vibrations  start  to  dominate  the  Fust  mode,  and  the  second  peak  shifts  upward 
and  it  may  eventually  move  out  of  the  operational  speed  range.  Therefore,  for  a  large  k 
value  associated  with  stiff  shafts  and  bearings,  a  single  degree  of  freedom  torsional 
model  as  used  in  Chapter  II  should  be  sufficient.  Another  effect  of  a  large  k  is  that  the 
ntp  subharmonic  solutions  replace  some  of  the  harmonic  solutions  as  shown  in  Figures 
3.28b,c  and  3.29b,c. 
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0  1  2 
frequency 


Figure  3.28.  ygia  versus.Qh,  plots  for  Fm  =0  1 .  Fah=0.05,  Fbj=0.  bh=bc=0.1. 

bbi=0.  a)  k*l,  b)  k=5,  c)  kalO;  n  harmonic,  A  period-one  non¬ 
harmonic,  o  period-n  snbharmonic,  •  quasi-periodic  or  chaos, - 

coiresponding  linear  response. 
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Figure  3.29.  ga  versus  Qh,ploisfor  Fm=0.1,  F.^0.05,  Fbi=0,  bh=bc=0.1,  bb,=0,  a) 
k  =  l,  b)  ka5,  c)  k=10;  o  harmonic,  A  period-one  non-harmonic,  o 

period-n  subharmonic,  •  quasi-periodic  or  chaos, - corresponding 

linear  response. 
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Figure  3.30.  y„ia  versus  £2h  plots  for  Fm=0.1,  Fbi?0,  bh=bc=0.1.  bbi=0  £33=0.05. 

Ci3=;23=0  0125.  ;,j=;22  =0.01,  k=l;  a)  Fah=0.05,  b)  Fah=1.0; 

harmonic,  o  period-n  subharmonic,  •  quasi-periodic  or  chaos, - 

corresponding  linear  response. 
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Figure  3.3 1 .  pa  versus^  plots  for  Fm=0. 1 ,  F^O.05,  Fbj=0,  bh=bc=0, 1 ,  1^=0,  a) 
k=l,  b)  k=5,  c)  k=10;  °  harmonic,  o  period-n  subharmonic,  •  quasi- 
periodic  or  chaos, - corresponding  linear  response. 
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Figures  3.30a,  3.31a  and  3.30b,  3.31b  examine  the  effect  of  load  ratio  for 
Fat>^m=0-^  and  1  respectively  with  Fm=0.1  and  k  =  1 .  The  frequency  response  is 
well-defined  and  dominated  by  the  harmonic  solutions  for  Ffth/Fm=0.5.  However 
when  Fflh/F,,,  is  increased  to  1,  the  region  beyond  0^=1. 0  becomes  strongly  non-linear 
consisting  of  subharmonic,  quasi-periodic  or  chaotic  responses.  This  was  also  seen  for 
the  gear  pair  model. 

B.  No  Gear  Backlash.  Non-linear  Bearings:  Now  consider  the  same  system  with 
bh=0  and  non-linear  bearings  described  by  the  approximate  model  given  in  Sections 
3.2.4.  Results  arc  presented  for  both  roller  and  ball  bearings  separately  since  the  radial 

clearance  in  a  ball  bearing  is  much  smaller  than  that  typically  found  in  a  roller  bearing 
with  the  same  inner  and  outer  race  diameters,  say  bbjro)(tr  *  10bbiM) .  For  each  case. 

bfcj  is  used  as  the  characteristic  length  bc. 

Figures  3.32  and  3.33  show  the  steady  state  frequency  response  plots  for  the 
roller  bearings  of  bj,j-0.01  mm  (i=l,2),  for  three  different  k  «=  kb,  / kh  values.  As 

shown  in  Figures  3.32a  and  3.33a  for  soft  bearings  (k  »=  1),  two  modes  of  interest  are 
weakly  coupled,  and  therefore  the  jump  phenomenon  is  seen  only  at  the  first  peak 
which  is  dictated  by  mostly  transverse  vibrations  Here,  the  bearing  non-linearities 
affect  the  transpose  JLplat  c.rt«*nts  y^j  «nd  y^  w^»ch  makes  the  first  mode  non-linear. 
The  second  peak  is  more  dependent  on  the  modal  interactions.  Here,  again  all  of  the 
solutions  are  purely  harmonic  and  only  no  impact  and  single  sided  impact  regimes 
exist.  However,  an  increase  in  k  introduces  chaotic  and  subharmonic  responses  as 
shown  in  Figures  3.32b,c  and  3.33b,c.  For  instance,  all  the  solutions  within  the  range 
0,6<f2h<1.7  are  non-harmonic  for  k  =  10;  the  frequency  response  is  no  longer  well 
defined  and  modal  coupling  is  sufficiently  strong  so  that  jump  discontinuities  are 
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observed  at  both  peaks. 

A  decrease  in  Fm/F^  enhances  the  degree  of  the  non-linearity  as  shown  in 
Figures  3.34a,b  and  3.35a,b.  At  Fm/F^sl  in  Figures  3.34b  and  3.35b,  the  responses 
within  0.5<ft|,<1.4  are  scattered  and  non-harmonic.  However  radial  preloads  Fbi  on 
bearings  should  help  the  mean  load  Fm  in  limiting  the  effect  of  non-linearities.  To 
prove  this  claim,  we  apply  a  high  radial  preload  to  both  bearings,  say  Fm  =  F^ ,  i  =  1,2 , 
corresponding  to  Fbj=0.25  for  Fm=1.0  in  the  dimensionless  form,  and  compare  results 
with  the  previous  case  of  Fbj-0  in  Figures  3.34a,c  and  3.35a,c.  Figures  3.34c  and 
3.26c  show  that  most  of  the  subharmonic  and  chaotic  solutions  of  Figures  3.34a  and 
3.35a  are  replaced  by  harmonic  solutions.  Hence  a  well  defined  frequency  response 
curve  with  clear  jump  discontinuities  is  obtained  by  applying  Fbi. 

Now.  replace  roller  bearings  with  ball  bearings  with  bc“bbj=0  .001  mm  and  again 

a 

investigate  the  effect  of  k  =  ky  /kh  in  Figures  3.36  and  3.37.  With  soft  bearings 

A 

(k  ss  1)  the  frequency  response  is  linear  (no  jumps)  which  indicates  that  there  are  no 
impacts  within  the  frequency  range  consideted  as  shown  in  Figures  3.36a  and  3.37a. 
However  at  a  larger  value  of  k ,  double  sided  impact  solutions  appear  since  the 
clearance  2bbj  is  very  small.  In  Figures  3.36b  and  3.37b  for  k  =  5,  the  frequency 
region  o?  0.65<f2b<0.8  consists  of  the  double  sided  solutions.  When  k  is  increased 
to  10,  the  same  behavior  is  also  seen  at  the  second  peak  as  shown  in  Figure  3.36c  and 
3.37c.  The  earlier  discussion  on  the  effect  of  Fjjh/Fm  and  Fbj/Fm  for  roller  bearings  is 
also  applicable  to  the  ball  bearings  as  well. 

3.6.  EXPERIMENTAL  VALIDATION 

Munro's  experimental  dynamic  transmission  error  results  [27],  which  were 
acquired  in  1962  using  a  four-square  spur  gear  test  rig,  have  been  used  earlier  to 
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frequency 


Figure  3.32.  y gja  versus  Qj,  plpts  for  the  case  of  roller^bearings;  Fm=1.0.  F^sO.S, 
^bi'O,  ^h— 0,  bjjj=bc=0.01mm,  a)  k=l,  b)  k=5,  c)  k=10;  c  harmonic,  A 
period-one  nonharmonic,  o  period-n  subharmonic,  •  quasi-periodic  or 
chaos,  corresponding  linear  response. 
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Figure  3.33.  pa  versus  plots  for  the  case  of  roller  bearings;  ^ah=0 -5. 

Fbi=0,  ^=0,  bb,=bc=0.01mm,  a)  k=l,  b)  k=5,  c)  k  =  10,  <=  harmonic,  A 
period-one  nonharmonic,  o  period-n  subharmonic,  •  quasi-periodic  or 
chaos, - corresponding  linear  response. 


frequency 

ygia  versus  Qh  plots  for  the  case  of  roller  bearings;  Fm=1.0,  1^=0, 
bbi=bc=0.01mm,  k=5;  a)  Fah=0.5,  Fbi=0,  b)  Fah=0.5.  Fbj=0.25,  c) 
Fah=F0,  Fbl=0;  <=>  harmonic,  A  period-one  nonharmonic,  o  period-n 

subharmonic,  •  quasi-periodic  or  chaos, - corresponding  linear 

response. 


frequency 

Figure  3.35.  pa  versus  Qh  plots  for  the  case  of  roller  bearings;  Fm  =  1.0,  bh=0. 

bbi=bc=0.01mm,  k=5;  a)  Fah=0.5,  Fbi=0,  b)  Fah=0.5,  Fbj=0.25,  c) 
Fah=^  ^>  °  harmonic,  A  period-one  nonhamionic,  o  period-n 

subharmonic,  •  quasi-periodic  or  chaos, - corresponding  linear 

response. 
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Figure  3.36. 


frequency 


yg,a  versus  £2b  plots  for  the  case  of  deep-groove^  ball  bearings; 
Fm  =  10,  Faj1=5,  Fbi=0,  bh=0,  bbj=bc=0. 001mm,  a)k=l,  b)k=5,  c) 
k=10;  “  harmonic,  A  period-one  nonharmonic  o  period-n  subharmonic, 
•  quasi-periodic  or  chaos, - corresponding  linear  response. 
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Figure  3.37.  p„  versus  12tl  plots  for  the  case  of  cieep-grooye  ball  bearings;  Fro=10, 
l‘ths3,  f;bj=0,  bh=0,  bf,j=bc=0.(K)l mm,  a)  k  =  ],  b)  k=5,  c)  k  =  10, 
harmonic,  A  period-one  nonharmonic,  o  period-n  subhatmonic,  •  quasi- 
periodic  or  chaos, - corresponding  linear  response. 


validate  the  spur  gear  pair  model  of  Chapter  II  assuming  that  the  modes  are  weakly 
coupled  with  each  other.  Now,  we  use  the  three  degree  of  freedom  non-linear  model  of 
Figure  3.2a  to  analyze  his  test  set-up  [27],  In  our  analysis,  bearings  are  assumed  to  be 
linear  since  they  were  highly  preloaded,  only  gear  backlash  non-linearity  is  considered 
with  2bjj=0.12  mm  [74].  Table  3.2  lists  the  system  parameters  used  in  equation  (3.3a) 
under  four  different  mean  loads.  The  system  damping  ratios  are  assumed  to  be  uniform 
at  each  load.  Average  mesh  stiffness  kh  and  alternating  load  ^ah  values  associated  with 
each  mean  load  Fm  ant  also  tabulated  in  Table  3.2. 

Figures  3.38  through  3.42  compare  dynamic  transmission  error  predictions  with 
the  measurements  [27].  At  the  design  load  (DL),  which  corresponds  to  the  minimum 
static  transmission  error,  good  agreement  is  seen  except  for  the  jump  discontinuity 
found  experimentally  at  the  second  peak.  In  our  model,  we  have  increased  slightly 
beyond  the  value  given  in  Reference  [27],  so  that  the  predicted  amplitude  of  the  first 
linear  peak  matches  very  well  with  the  experiment.  Such  slight  changes  in  critical 
system  parameters  such  as  F^  and  although  they  are  varied  within  the  experimental 
uncertainties,  may  alter  the  frequency  response  drastically  as  illustrated  in  Figures 
3.39a  and  3.39b.  In  Figure  3.39a  we  note  that  when  we  reduce  the  force  ratio  Fm/F^ 
from  30  (the  value  given  in  Reference  [27])  to  10,  a  jump  discontinuity  at  the  second 
peak  is  seen  which  is  compatible  with  the  experiment.  Similarly,  a  small  change  in  the 
damping  values,  which  are  not  reported  in  Reference  [27],  also  affects  the  frequency 
response  as  shown  in  Figure  3.39b.  Figures  3.40  to  3.42  compare  the  dynamic 
transmission  error  spectra  at  3/4,  1/2  and  1/4  of  the  design  load.  From  these  spectra  we 
conclude  that  our  proposed  theory  agrees  well  with  the  experimental  results  of  Munro 
[27]  bo. i  inalitatively  and  quantitatively,  although  amplitudes  in  the  off-resonance 
regions  are  slightly  off.  Such  amplitudes  are  close  to  the  noise  floor  in  the  experiment. 
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Table  3.2.  Parameters  of  the  Munro's  experimental  set-up  extracted  from  Reference 


[27], 


Mean  load,  Fm 

design  load(DL) 

0.183 

3/4  of  DL 

0.146 

1/2  of  DL 

0.105 

1/4  of  DL 

0.0579 

Fah 

0.0058t 

0.0178 

0.0296 

0.0393 

kb  (N/m) 

1.16x109 

U6xl09 

1.16x109 

1.16x109 

kh(NAn) 

3.44x108 

3.22x108 

3.01x10s 

2.72x10s 

KU-  k22 

0.950 

0.966 

0.983 

1.007 

k13*  *23 

0.242 

0.242 

0.242 

0.242 

Cll-  ^22tf 

0.01 

0.01 

0.01 

0.01 

Cl 3  •  C23tf 

0.00375 

0.00375 

0.00375 

0.00375 

C33n 

0.015 

0.015 

0.015 

0.015 

t  modified  so  that  the  linear  peak  matches  with  the  measured  value, 
tt  estimated  using  the  experimental  data  of  Reference  [27]. 
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Figure  3.38.  Comparison  of  theory  with  Munro's  [27 ]  experiment  at  the  design  load. 


158 


a) 


o 

t 


b) 


frequency 


Figure  3.39.  Effect  of  a)  alternating  load  and  b)  damping  on  the  dynamic  transmission 
error  spectrum  at  the  design  load. 
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Further,  discrepancies  between  the  theory  and  the  experiment  may  be  due  to  the 
assumptions  made  in  this  study.  It  should  however  be  pointed  out  that  the  experiment, 
even  though  it  was  very  precise,  was  conducted  almost  three  decades  back  before  the 
advent  of  modem  non-linear  dynamics  and  chaos  science  [66-71],  Interestingly, 
Munro  [27]  had  also  reported  subharmonic  responses  at  3/4,  1/2  and  1/4  of  the  design 
load,  and  "unrepeatable  responses"  (probably  chaos)  at  1/4  of  the  design  load.  Our 
results  agree  with  such  experimental  observations  as  well.  Both  experiment  and  theory 
show  that  the  chaotic  solutions  dictate  the  frequency  response  beyond  the  second 
resonance  for  1/4  of  the  design  load  in  Figure  3.42. 

3.7.  CONCLUSION 

This  analytical  study  on  the  non-linear  dynamics  of  a  geared  rotor-bearing  system 
with  gear  backlash  and  bearing  clearances,  as  excited  by  the  internal  static  transmission 
error  and/or  external  torque  pulsations,  has  made  a  number  of  contributions  to  the  state 
of  the  art.  First,  an  approximate  non-linear  model  of  the  rolling  element  bearings  with 
clearances  has  been  developed  and  validated.  Second,  our  mathematical  model  has 
been  validated  as  it  compares  well  with  a  previous  experimental  study,  and  several  key 
parameters  such  as  the  mean  load,  mean  to  alternating  force  ratio,  radial  bearing 
preload,  bearing  stiffness  and  damping  have  been  identified.  Third,  conditions  for  the 
chaotic  and  subharmonic  resonances,  and  the  routes  to  chaos  have  been  identified. 
Forth,  non-linear  modal  interactions  have  been  examined.  Finally,  on  a  more 
fundamental  note,  our  study  enriches  the  current  literature  on  coupled  vibro-impact 
pairs  [52,54,55i. 
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CHAPTER  IV 


INTERACTIONS  BETWEEN  TIME-VARYING  MESH  STIFFNESS 
AND  CLEARANCE  NON-LINEARITIES 

4.1.  INTRODUCTION 

Dynamic  models  of  geared  systems  can  be  classified  into  four  main  groups.  The 
first  group  includes  ''rear  time-invariant  (LTI)  models,  as  evident  from  an  extensive 
review  of  the  literature  given  in  Reference  [14].  The  second  group  considers  linear 
time-varying  (LTV)  mesh  stiffness  kj,(t)  in  the  analysis  [27,29.62].  Periodic  variation 
in  kh(0  is  due  to  the  changing  of  the  number  of  conjugate  teeth  pairs  in  contact  during 
the  convolute  action.  Accordingly,  the  system  is  excited  parametrically  as  well  as  by 
the  static  transmission  error  e(t)  introduced  by  kinematic  errors  and  tooth  deflections. 
In  this  case,  the  equation  of  motion  of  the  gear  pair  essentially  reduces  to  the  Mathieu's 
or  Hill's  equation  with  a  periodic  external  forcing  function.  The  third  group  includes 
gear  backlash  in  the  models,  but  with  time-invariant  average  mesh  stiffness  k^  *  kh(t) 

[25,30,31,33,34,41,42].  It  should  be  noted  that  backlash  is  bound  to  exist  either  by 
design  or  due  to  manufacturing  errors  and/or  wear  in  any  gear  pair.  Finally,  the  last 
group  considers  be  h  gear  backlash  and  mesh  stiffness  variation  simultaneously 
[24,35-37,53].  However,  none  of  these  studies  have  addressed  explicitly  the  effect  of 
kh(t),  including  its  interaction  with  the  backlash  non-linearity,  on  the  steady  state 
frequency  response.  Ozguven  and  Houser  [24]  have  attempted  to  analyze  this  problem 
by  replacing  kh(t)  with  a  constant  mesh  stiffness  and  by  defining  the  "loaded  static- 
transmission  error"  excitation  at  the  mesh  point.  But,  Ozguven  [75]  has  stated  recently 
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that  this  approach  may  not  work  depending  on  the  system  parameters,  and 
recommended  a  detailed  investigation  of  this  issue. 

Numerous  publications  have  analyzed  LTV  systems  [76-78]  and  time-varying 
systems  with  quadratic  and/or  cubic  non-linearities  [79-82].  However,  such  studies  are 
not  directly  applicable  to  the  geared  rotor  systems  which  are  the  main  focus  of  this 
chapter.  This  problem  requires  the  solution  to  a  set  of  time-vary  ing  differential  equations 
with  clearance-type  non-linearities  as  excited  by  a  periodic  force  generated  at  the  gear 
mesh. 

4.2.  PROBLEM  FORMULATION 

A  generic  geared  rotor-bearing  system,  which  consists  of  a  spur  gear  pair 
mounted  on  flexible  shafts  which  are  supported  by  rolling  element  bearings  as  shown 
in  Figure  4.1a,  is  considered  here.  The  gear  box  is  assumed  to  be  rigid.  The  effect  of 
the  prime  mover  or  the  load  inertia  is  not  considered  assuming  that  such  inertial 
elements  are  connected  to  the  gear  box  through  soft  torsional  couplings.  Further,  we 
assume  that  the  system  is  symmetric  about  the  plane  of  the  gears  and  the  axial  motion 
(parallel  to  the  shafts)  is  negligible.  The  governing  equations  of  motion  can  be  given  in 
the  matrix  form  as 

[M]{q"(t)}  +  [c]{q'(t)}  +  [K(t)]{f(q(t))}  =  {F(t )}  (4.1) 

where  [m]  is  the  time-invariant  mass  matrix  and  {q(t)}  is  the  displacement  vector. 
Here,  damping  matrix  [c]  is  assumed  to  be  time-invariant,  as  the  effect  of  the  tooth 

separation  and  time-varying  mesh  properties  on  mesh  damping  are  considered 
negligible;  validity  of  this  assumption  will  be  examined  later.  The  stiffness  matrix 
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a)  Generic  geared  rotor-bearing  system,  b)  clearance  non-linearity  in 
gears  and  bearings. 
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[K(t )]  is  considered  to  be  a  periodically  time-varying  matrix  given  by 
[K(t)]  =  [K(t  +  2rc/Qh)j  where  Qh  is  the  fundamental  gear  mesh  frequency  The 
non-linear  displacement  vector  {f (q(t ))}  includes  the  radial  clearances  in  bearings  and 
the  gear  backlash  as  shown  in  Figure  4.1b,  and  the  forcing  vector  {F(t )}  consists  of 

both  external  torque  and  internal  static  transmission  error  excitations. 

This  chapter  extends  our  previous  non-linear  single  degree  of  freedom  spur  gear 
pair  model  of  Chapter  II  and  multi-degree  of  freedom  geared  rotor-bearing  system 
model  of  Chapter  ID  by  including  time-varying  mesh  stiffness  kh(t),  and  investigates 
its  effect  on  the  frequency  response  of  lightly  and  heavily  loaded  geared  systems. 
Interactions  between  mesn  stiffness  variation  and  system  non-linearities  associated  with 
gear  backlash  and  radial  clearances  in  rolling  element  bearings  are  also  considered. 
Resonances  of  the  corresponding  LTV  system  associated  with  the  parametric  and 
external  excitations  are  identified  using  the  method  of  multiple  scales.  Our  formulation 
will  be  validated  by  comparing  predictions  with  available  experimental  results  [27,28] 

4.3.  MATHEMATICAL  MODEL 

A  reduced  order  form  of  the  multi-degree  of  freedom  system  of  equation  (4.1), 
the  three  degree  of  freedom  non-linear  model  of  the  geared  rotor  system,  which  has 
been  used  in  Chapter  III,  is  considered  here.  This  includes  gear  inertias  Igl  and  Ig2, 
gear  masses  mgj  and  and  base  circle  diameters  d-i  and  dg2,  as  shown  in  Figure 
4.2a.  The  gear  mesh  is  described  by  a  non-linear  displacement  function  fh  with  time 
varying  stiffness  kh(t)  and  linear  viscous  damping  c^  Bearings  and  the  shafts  that 
support  the  gears  are  modeled  by  equivalent  elements  with  viscous  damping 
coefficients  c^i  and  C(,2  and  non-linear  springs  defined  by  force-displacement  functions 
fbl  and  f{,2  which  are  approximated  by  piece-wise  linear,  dead  zone  type  non-linearities 
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as  suggested  in  Chapter  HI.  Both  low  frequency  external  excitation  due  to  torque 
fluctuations  and  high  frequency  internal  excitation  due  to  the  static  transmission  error 
e(t)  are  considered  in  the  formulation.  Input  torque  fluctuation  is  included,  but  the 
output  torque  is  assumed  to  be  constant,  i.e.  Tgi(t)  =  Tgjm  +  Tgja(t)  and 
Tg2(t)  =  Tg2m.  External  radial  prek  acs  and  are  also  applied  to  both  rolling 
element  bearings. 


4.3.1.  Equations  of  Motion 

Equations  of  coupled  transverse-torsional  motion  of  the  geared  rotor-bearing 
system  shown  in  Figure  4.2a  with  the  displacement  vector 
{q(t)}  =  |ygi(t),  yg2(t),  p(t)j  are  given  in  matrix  form  as  follows: 
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Figure  4  2  Nonlinear  models  of  the  yea;ed  system  of  Figure  !;  a)  three  degree  of 
freedom  system,  b)  gear  pair  model  with  a  single  degree  of  freedom. 
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Here,  (  )'  donates  derivative  with  respect  to  time  t ,  ygj  and  0gi  are  the  transverse  and 
torsional  displacements  of  the  i-th  gear  (i=l  ,2),  mcj  is  the  equivalent  gear  pair  mass, 
Fm  is  the  average  force  transmitted  through  the  gear  mesh,  and  Fay(t)  is  the  fluctuating 
force  related  to  the  external  input  torque  excitation.  Equations  (4.2a-h)  have  been  used 
previously  in  Chapter  III  except  for  kb(t )  term,  which  is  expanded  in  equations  (4.2a- 
h)  in  the  Fourier  series  form.  Here  p(t)  is  the  difference  between  the  dynamic 
transmission  error  and  the  static  transmission  error  e(t).  A  dimensionless  form  c. 
equation  (4.2a)  is  obtained  by  letting  ygj(0  =  yg/t)/bc ,  p(t)  =  p( t ) /  bc, 
0)n  =  -v/khn,  /  mc| ,  tobl  =  y'kbj  /  .tigj  (i=l,2)  and  t  =  (0„t  where  bc  is  the  characteristic 
length.  Here,  we  consider  periodic  excitation  for  both  e(t)  and  FaT(t)  as 
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c(t)=  J-£rcos  (rft(,i  +  4>Cr)-  ^al(‘)=  I^aTr  cos^rftyt  +  4fj,  ]  where  ftb  and  ftr-  are 
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the  fundamental  excitation  frequencies  of  internal  displacement  and  external  torque 
fluctuations,  respectively.  Further,  define  dimensionless  excitation  frequencies 
ft|,  -  fth  /  ton  and  ftj=ftx/(un  to  yield  the  following  governing  equations  of 
motion  in  the  dimensionless  fonn; 
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where  Fj,j  (i=l,2)  and  Fm  are  the  dimensionless  components  of  the  mean  force  vector 
{F)m,  and  Fa-r(t)  and  F^t)  pertain  to  the  alternating  external  excitation  |F(t))e  and 
internal  excitation  ( F(t))j  force  vectors,  respectively. 

4.4.  GEAR  PAIR  STUDIES 

First,  we  consider  the  spur  gear  pair  model  of  Chapter  II,  as  shown  in  Figure 
4.2b,  and  investigate  the  effect  of  kh(t)  on  the  steady  state  frequency  response.  The 

equation  of  motion  of  the  gear  pair  is  obtained  from  equation  (4.3a)  by  substituting 
ygi(t)=0  for  fixed  gear  centers.  Neglecting  the  input  torque  variations  for  the  sake  of 
convenience,  Faj(t)=0,  we  obtain  the  following  equation  of  motion  of  a  loaded  gear 
pair  with  time-varying  mesh  stiffness  and  backlash,  as  excited  by  the  static 
transmission  error. 


p(t)  +  2£33p(t)  +  K33(t  )fh(p(0)  =  Fm  +  Fah(0 


(4,4) 


where  £33,  <33(0,  fh(p),  Fm  and  F^  are  still  given  by  equation  (4.3).  Since  non?;  c, 
the  existing  analytical  solution  methods  arc  found  to  be  suitable  for  this  problem,  the 
digital  simulation  technique  is  used  to  solve  equation  (4  4)  A  'jth-bth  order  Rungc 
Kutta  numerical  integration  algorithm  w;ti.  variable  ;me  step  |50J  is  used  hern.  This 
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technique  has  already  been  employed  successfully  in  Chapter  D  for  the  gear  pair 
problem  with  time-invariant  mesh  stiffness. 

First  consider  only  sinusoidally  varying  mesh  stiffness  and  the  static  transmission 
error  excitation,  i.e.  K33O)  =  1  +  Ej  cos(Qbt  +  4>bl  >  where  £=£j  and 
Fjjjft)  =  Fjtjjfttj2  cos(nbt  +  4»el).  The  tooth  deflection  or  F^ft)  is  maximum  under  the 
applied  mean  load  Fm  when  K33(t)  is  minimum  (i.e.  minimum  number  of  gear  pair  in 
contact).  Similarly,  minimum  F^ft)  corresponds  to  maximum  >€33(0.  Therefore,  there 
is  an  out  of  phase  relationship  between  F^t)  and  <33(1),  i.e.  4>bl=<t>el+,t-  Here*  we 
set  <|>b|=7c  and  4>ci=0  for  the  sake  of  convenience.  Hence,  equation  (4.4)  is  modified 
to: 


p(t)  +  2^33p(t)  +  [l-£cos(nht)]fh(p(t))  =  Fm  +  FabQbcos(ftht).  (4.5) 

4.4.1.  Linear  Time-varying  (LTV)  System 

For  zero  gear  backlash  bj,,  the  gear  mesh  displacement  function  is  fb(p)=p  Hence 
equation  (4.5)  reduces  to  an  LTV  equation  as  follows: 

p(t)  +  2;33p(t)  +  [l-£cos(f2ht)]p(t)=  Fm  +  Fabnbcos(ftbt)  (4.6) 

Approximate  analytical  solutions  of  similar  LTV  differential  equations  with  parametric 
and  external  excitations  are  already  given  in  the  literature  [76.77.79-82J.  Therefore,  we 
will  not  attempt  to  solve  equation  (4.6)  completely.  Instead  we  identify  the 
corresponding  resonances  using  the  method  of  multiple  scales  (76).  A  first  order 
uniform  solution  is  given  by  an  expansion  having  the  following  form  where  the  scalar 
parameter  e  is  assumed  to  be  very  small. 
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(4.7) 


p(t,e)  =  p0(T0,Ti )  +  £p1(T0,Tj)  +  0(e2  ) 

where  TD  =  ent  is  the  n-th  time  scale.  Expansions  for  the  derivatives  with  respect  to  t 
is  obtained  in  terms  of  the  partial  derivatives  with  respect  to  time  scales  Dn  =  d  /  dTn 

d  /  dt  =  D0  +  el>!  +  0(£2 ) ;  d2/dt2  =  D2+2eD0D,+0(E2)  (4.8a, b) 

Substituting  equations  (7)  and  (8)  into  equation  (6)  and  equating  like  powers  of  e,  with 
the  external  force  being  applied  at  0(E°),  one  obtains 

D0Po+Po  =  Fm+Fahlahcos(QhTo)  (4.9a) 

d2pi  +  Pi  =-2noD,p0-pD0p0  +  p0cos(flhT0)  (4.9b) 

where  £^2^33.  The  solution  of  equation  (9a)  is  given  in  the  complex  domain  as 

Po  =A(T,)eiT°  +  Fm  +  Acii2hT°  +  cc,  A  =  /  2(1  -ftg)  (4.10a.b) 

where  cc  represents  the  complex  conjugate  terms  and  i  =  V-T .  At I2h  « 1 ,  we  note  that 
(l-ftfl,  will  make  the  amplitude  of  the  response  boundless;  this  is  the  primary 
resonance  since  it  appears  in  the  first  order.  Substituting  equation  (4.10)  into  equation 
(4.9b)  we  get 
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DoPl  +P1  =  -2iD1AelT°  -i^Ae,T°  -ipAnhe'nhT°  +  ^  Ael(1+nh)T° 

+  -Acl(l-tlh)T°  +-Ae2,nhT°  +  —  Froe^1'1®  +  A  +  cc  (4.11) 
2  2  2  m 

Besides  the  primary  resonance  at  Dtj,  "■  1 ,  the  particular  solution  of  equation  (4.11)  ha3 
other  secular  terms  when  2fth  *  1  and  *>  2 .  At  20),  •  1,  summation  of  the 
external  excitation  frequency  and  the  parametric  excitation  frequency  is  close  to  the 
dimensionless  natural  frequency  which  is  unity.  One  could  also  observe  resonances  at 
Ofl  *n,  n>2  when  higher  scales  are  considered  in  equation  (4.7).  Jn  nummary, 
equation  (4.6)  has  resonances  at  0|,  *  0,5  and  n,  where  n=l  ,2.3... 

Now,  we  solve  equation  (4.6)  using  digital  simulation.  Figure  4.3  shows  the 
steady  state  frequency  response  curves  pa(Oh)  and  pm(Oj,)  for  a  lightly  loaded  system 
with  Fm=0.1 ,  FahJ=0.05,  £33-0-05,  b),=0  and  four  different  e  values.  Note  that  e-0 
represents  the  LTI  system.  For  e>0.  we  do  not  observe  any  multi-valued  regions  and 
jump  phenomenon,  similar  to  those  seen  for  the  non-linear  systems.  Mesh  stiffness 
variation  e  has  a  negligible  effect  on  the  natural  frequency  which  corresponds  to  the 
largest  peak  in  Figure  4.3a;  but  an  increase  C  in  amplifies  p„  in  both  resonance  and  off- 
resonance  regions.  Here  pm  is  no  longer  uncoupled  from  pa,  and  it  varies  substantially 
in  the  vicinity  of  primary  resonance  at  flj,=  l,0,  and  the  parametric  resonance  at 
fth=0.5,  especially  for  a  large  e  confirming  the  analysis  given  earlier.  Note  that  our 
study  is  limited  to  only  these  two  resonances  as  an  investigation  of  subharmonic 
resonances  at  ft;,  *  n,  n>l  is  beyond  the  scope  of  this  study. 
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Figure  4.3.  Frequency  response  spectra  of  a  lightly  loaded  LTV  gear  pair  with 
sinusoidal  e(t)  and  khO).  Fm=0.1,  F^cO-OS,  £$3=0.05,  bh=0  and  four 
different  e  values;  a)pa  versus  Hjj,  b)  pm  versus  fljj. 
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4.4.2.  Non-linear  Time-varying  System 

Next  we  consider  the  time-varying  non-linear  system  with  backlash,  given  by 
equation  (4.5).  Figure  4.4  shows  pa(£2b)  and  pm(£2b)  spectra  with  Fm=0.1,  FaJ=0.05, 
^33=0.05,  and  four  £  values.  In  this  case,  we  notice  a  jump  discontinuity  at  the 

resonant  frequency  and  a  dual-valued  region  bounded  by  jump-up  and  jump-down 
transition  frequencies.  For  e=0,  pm  *pm(flb)  is  evident  in  the  no-impact  regime, 

similar  to  the  results  given  in  Chapter  II.  But  when  E>0,  the  transition  frequencies 
which  define  the  jump  phenomenon  become  smaller,  and  pa  at  the  jump-down 
frequency  grows  with  increasing  e.  This  indicates  that  kb(t)  enhances  the  degree  of 

non-linearity  associated  with  the  gear  backlash.  Similar  to  the  results  of  LTV  system  in 
Figure  4.3,  we  again  observe  in  Figure  4.4a  a  parametric  resonance  at  £^=0  5,  which 
is  strong  for  £=0.2  and  0.4  curves.  Figure  4.5  compares  the  time  histories  p(t)  for 
each  e  values  at  the  parametric  resonance  peak,  Qb3^  For  LTI  system  (£=0),  a 
harmonic  no-impact  type  steady  state  solution  exists  as  illustrated  in  Figure  4.5a.  With 
increasing  E,  this  solution  is  transformed  into  a  non-harmonic  periodic  solution  with  a 
larger  peak  to  peak  value  as  a  result  of  the  parametric  mesh  stiffness  excitation,  as 
shown  in  Figures  4.5b-d. 

Now,  consider  a  heavily  loaded  gear  pair  with  Fm=0.1,  Faj=0.01  and  £33=0.05. 
Figures  4.6a  and  4.6b  show  pa(£2h)  and  pm(£2h)  spectra,  respectively.  In  this  case, 
tooth  separation  does  not  occur  for  £=0  which  results  in  a  linear  frequency  response 
curve  with  a  constant  pm.  A  small  jump  is  seen  at  e=0. 1 ,  and  this  jump  becomes  larger 
for  £=0.2  and  0.4;  overall  alternating  amplitudes  pa  at  resonance  and  within  the  off- 
resonance  regions  increase  considerably.  This  suggests  that  the  mesh  stiffness 
variation  is  especially  important  for  a  heavily  loaded  gear  pair  with  backlash. 
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Figure  4.4.  Frequency  response  spectra  of  a  lightly  loaded  non-linear  gear  pair  of 
Figure  2bwith  sinusoidal  e(t)  and  kb(t),  Fm-0.1,  FjjjjsO.05,  £33=0.05, 
bij=bc  and  four  different  e  values;  a)pa  versus  b)  pm  versus 
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Figure  4.6.  Frequency  response  spectra  of  a  heavily  loaded  non-linear  gear  pair  with 
sinusoidal  e(t)  and  kh^O,  Fm=0.1,  Fjhj^.01,  £33=0.05,  bh=bc  and  four 
different  e  values;  a)pa  versus  b)  pm  versus  £2^. 
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4.4.3.  Non-linear  Time-varying  Mesh  Damping 

In  order  to  investigate  the  effect  of  non-linear,  time-varying  gear  mesh  damping 
on  the  steady  state  frequency  response,  we  consider  a  sinusoidally  varying  gear  mesh 
damping  £33(0  which  is  assumed  to  be  proportional  to  the  mesh  stiffness  function 
*33(t),  and  a  non-linear  velocity  function  gh  associated  with  tooth  separation.  Hence 
the  governing  equation  of  motion  is 

p(t )  +  2^33(t)gh(p(t))  +  [1  -  Ecos(Qht)]fh(p(t))  =  Fm  +  Fahftb  cos(Qht);  (4. 12a) 


^33^)"  C33k33(*)  =  £33^  +  £COS(Qht  +  <|>hi))  ; 


(4.12b) 


gh(p('))  = 


0; 

p(t); 


bh  b|i 

—  <  p  <  — 

bc  K  bc 

else 


(4.12c) 


Figure  4.7  compares  frequency  response  spectra  for  the  cases  of  a)  non-linear  time- 
varying  gear  mesh  damping  defined  by  equation  (4.12),  and  b)  linear  time-invariant 
mesh  damping,  £33  *£33(0  with  gh(p(0)  =  p(t),  as  given  by  equation  (4.5).  As 

shown  in  Figure  4.7,  both  mesh  damping  models  yield  virtually  the  same  spectra. 
Therefore  a  linear  time-invariant  gear  mesh  damping  model  can  be  used  without  losing 
any  accuracy. 

4.4.4.  Periodic  e(t)  and  kj,(t) 

Up  to  now,  we  have  only  considered  the  sinusodally  varying  kh(t )  and  e(t ).  In 
real  geared  systems,  however,  both  kh(t)  and  e(t)  are  periodic  which  can  be 
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Figure  4.7.  Frequency  response  spectra  of  gear  pair  corresponding  to  non-linear  time 
varying  and  linear  viscous  damping  models;  a)pa  versus  b)  pm  versus 
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expressed  In  ihc  Fourier  series  form  Figure  4  8  illustrates  k,,(t)  and  e(t)  for  a  low 
contact  ratio  Hpiii  gem  pin,  thrhc  o!  lime  histories  arc  predicted  using  an  existing 
spur  grttt  tooth  model  |2l  ,  'nr**  ftm  six  Fourier  coefficients  of  kh(i)  and  e(i)  arc 
given  In  Table  4,1,  Note  that  ii-ili  harmonic  of  l(i)  causes  a  resonance  around 
12,,*  I/m.  Hence  higher  harmonics  are  important  only  at  low  frequencies,  Accordingly, 
only  the  first  three  Fourier  coefficient*  of  e(t)  are  deemed  sufficient  h«:re,  given  the 
naluie  of  the  practical  gear  nolle  and  vibration  problems,  Similarly.  k,,(t)  is  also 
liuiiutlecJ  foi  the  sake  ol  convenience,  retaining  the  first  three  terms,  Higher  terms 
(n>3)  could  l»c  easily  included  in  the  analysts  if  ntcexsaiy.  Therefore  the  governing 
equation  with  I  rt|<u-o  is 


pirn  2£'«M)p(  !H  |!  -  Vt|  ios(!l2,,uJf|,(p(IH  =  F,„  ♦  VF,flfni2„rcos(ii2,,i) 

M  i  -I 


(4.13) 


Figutes  4 Vh  and  4  9 h  show  ficquency  icsponsc  spectra  foi  Fm*0.1,  Fi|lj"0.()3. 
f,I!,2»0  02,  Fi|(c*U  01,  i  j/rj*2,  r,/t  »*4  and  fom  different  c , «*£  values 

As  shown  in  Figtti?  4,9,  two  oilier  jump  discontinuities,  which  are  clearly  associated 
with  1^,2  and  F,j,\,  arc  found.  These  Jumps  me  similar  to  the  fundamental  resonance 
|*eak  except  they  have  lower  amplitude*  In  order  for  these  two  Jumps  to  exist, 
and  I'nj,  i  excitation*  have  to  l*e  large  Once  again,  n  i*  evident  that  the  periodic  mesh 
allflness  enhances  the  extent  of  backlash  non-liueutily 
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Figure  4.8  e(t)  and  ^^(1)  plots  for  a  low  contact  ratio  spur  gear  pair  [21].  Here, 
number  of  teeth  is  28,  diametral  pitch  is  8,  the  pressure  angle  is  20 
degrees,  and  a  6x1  O'4  in.  tip  modification  starting  at  26  degrees  roll 
angle  is  applied. 


Tabic  4.1  Fourier  coefficients  of  e(t)  and  given  in  Figure  4.8 


Harmonic  n 

e(t) 

4im) 

Amplitude 

M> 

(MN/m) 

0 

20.00 

107.8 

1 

3.45 

10.3 

2 

0.26 

5.7 

3 

0.45 

1.3 

4 

0.38 

4.0 

5 

0.27 

1.5 

6 

0.26 
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Figuie  4,9,  Frequency  response  spectra  of  a  non-linear  gear  pair  with  periodic  e(t )  and 
kh(t).  Fm-0.1,  F.|, 1-0.03,  Flh2-0.02,  F^-O-Ol.  $w*0.05,  bh*bc, 
t  j/t;2=2.  tj/Cj-4  and  four  different  Ej»i  values;  a)p„  versus  Oh,  b)  pni 
versus  il|, 
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4.5.  GEARED  ROTOR-BEARING  SYSTEM  STUDIES 


Next,  we  solve  the  three  degree  of  freedom  non-linear  model  of  Figure  4.2a  and 
equation  (4.3a)  using  the  digital  simulation  technique  for  excitation  Fah(t)  with 
FaT(t)=0. 


4.5.1.  Sinusoidal  e(t)  and  kh(t) 

Similar  to  Section  4.4,  assume  sinusoidal  forms:  <33(0  =  l+ecos(Qht  -Mt>h|) 
and  Fjhlt)  =  F^Of,2  cos(Oht  +  ^el  )•  Equation  (4.3a)  reduces  to  an  LTV  matrix 
equation  when  gear  backlash  and  radial  bearing  clearances  are  set  to  zero,  i.e. 
fhOW))5*^1)  and  fbi(ygi(,))=>'g«(,)>  i=1*2-  Figure  4.10  shows  pa(Dh)  30(1  >’gia(^h) 
spectra  of  this  LTV  system  under  heavily  loaded  conditions  with  Fm=0.1,  Fahl=0.0l, 
C>33=0. 05,  ^3=0.0125  and  ^=0.01,  i=l,  2.  Note  that  ygJa  =yg2a  since  the  gear  ratio 
is  one.  As  shown  in  Figure  4.10,  frequency  response  has  peaks  at  nh=coj=0.4  and 
Qh=0)jj=1.25  where  toj  and  a>u  are  the  natural  frequencies  of  the  LTI  system, 
corresponding  to  the  first  two  coupled  transverse-torsional  modes.  In  Figure  10a,  we 
observe  a  parametric  resonance  at  2L2^=(Oi  sitndar  to  the  gear  pair  model  of  Section  4.4 
whereas  the  parametric  resonance  at  2f2h=(Ojj  is  more  obvious  in  pa(Qh)  spectrum  of 
Figure  4.10b  for  a  larger  l  value.  Again  like  the  LTV  gear  pair  results  of  Figure  4.3, 
ygla  and  pa  are  amplified  with  increasing  E  in  both  resonance  and  off-resonance 
regions. 

Now  introduce  the  gear  backlash  bh=bc  to  this  heavily  loaded  system; 
corresponding  results  are  given  in  Figure  4.1 1.  The  most  significant  effect  of  kj,(t)  is 
that  it  interacts  with  the  gear  backlash  non-linearity  to  develop  a  jump  discontinuity  at 
the  second  primary  resonance  peak.  This  is  clearly  evident  from  the  fact  that  either 
backlash  or  k|,(t)  alone  can  not  cause  a  jump  for  such  a  heavily  loaded  system. 
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Figure  4.10.  Frequency  response  spectra  of  a  LTV  geared  rotor-bearing  system  of 
Figure  2a  with  sinusoidal  e(t)  and  kh(t),  Fm=0.1,  F^isO.Ol,  £33=0.05, 
£33=0.05,  £i3=0.0125  and  £jj=0.01,  k„=0.5,  bjj=bhbi=0,  i=l,  2,  and  four 
Cpe  values;  a)  ygja  versus  i2b,  b)  pa  versus  Qb. 
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Besides  an  increase  in  amplitudes  in  Figures  4.11a  and  4.1  lb,  chaotic  and/or  quasi- 
periodic  solutions,  which  do  not  exist  in  heavily  loaded  time-invariant  systems  as 
shown  in  Chapter  HI,  are  also  predicted  in  the  following  regimes:  1 .35<£ih<l  .65  for 
£=0.2  and  1.25<Qb<1.7  for  e=0.4.  Similar  to  the  LTV  system  of  Figure  4.10, 
parametric  resonances  at  2f2b=(Oj  and  2£2h=(Dn  are  again  observed  here. 

Next,  radial  clearances  in  bearings  bbj=bc  are  considered  for  a  heavily  loaded 
system  with  zero  gear  backlash,  Fm=1.0,  F^sO-l,  £33=0. 05,  £13=0. 0125  and 
£ii=0.01,  i=l,  2.  In  this  case,  mesh  stiffness  variation  and  non-linearity  exist  in 
different  components.  As  shown  in  Figure  4.12,  £  again  increases  the  amplitudes,  but 
the  spectrum  shape  is  essentially  the  same.  This  is  further  evident  from  Figure  4.13, 
given  for  a  lightly  loaded  drive  with  Fm=1.0  and  F^isO.5.  This  suggests  a  weak 
interaction  between  the  gear  mesh  stiffness  variation  and  bearing  non-linearities. 

4.5.2.  Periodic  e(t)  and  kb(t) 

Similar  to  the  gear  pair  analysis  of  Section  4.4.4,  we  consider  periodic  K33O)  and 
Fah(t)  with  three  Fourier  coefficients.  Figure  4.14  shows  spectra  for  Fm=0.1, 
Fahi=0  01,  Fab2=0.004,  Fah3=0.002,  C33=0.05,  ;i3=0.0125  and  <^=0. .01,  i=l,  2; 
£i/£2=2,  £i/£2=4,  ard  Ej=e=0,  0.1  and  0.2.  Here  only  the  gear  backlash  non-linearity 
(bb=bc  and  bbj=0)  is  considered.  Peaks  at  f2b=G>j/n  and  nb=tO]j/n,  n=l,2,3  are 
predicted  corresponding  to  excitation  F^.  Periodic  mesh  stiffness  enhances  the 
alternating  amplitudes  over  the  values  given  by  the  time-invariant  mesh  stiffness  case 
and  introduces  a  jump  discontinuity  at  the  second  primary  resonance  similar  to  the 
sinusoidal  mesh  stiffness  case. 
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Figure  4. 1 1 .  Frequency  response  spectra  of  a  non-linear  geared  rotor-bearing  system 
with  sinusoidal  e(t)  and  kb(t),  Fm=0.1,  Fahl=0.01,  C33  =  0.05, 
^33=0.05,  ^3=0.0125  and  ^=0.01,  ^=0.5,  bb=bc,  bt,j=0,  i=l ,  2,  and 
four  £j=E  values;  a)  y^  versus  b)  pa  versus  Hb. 
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Figure  4.12.  Frequency  response  spectra  of  a  non-linear  geared  rotor-bearing  system 
with  sinusoidal  e(t)  and  kh(*).  Fm=l .0,  F^isO-l,  £33=0.05,  £33=0.05, 
£i3=0.0125  and  £jj=0.01,  <^=0.5,  b^^bc,  bjj=0,  i=l,  2,  and  four  £j=e 
values;  a)  y^a  versus  Qj,,  b)  Pa  versus  flh. 
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Figure  4.13.  Frequency  response  spectra  of  a  non-linear  geared  rotoi -be aring  system 
with  sinusoidal  e(t)  and  kb(0.  Fm=1.0,  F^jsO.S,  £33=0.05,  ^33=0.05 , 
£i3=0.0125  and  ^-0.01 ,  ku=0.5,  bbi=bc,  ^=0,  i=l.  2,  and  four  e,=e 
values;  a)  yg,a  versus  £!{,•  P)  Pa  vcrsus 
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Figure  4.14.  Frequency  response  spectra  of  a  non-linear  geared  rotor-beaiing  system 
with  periodic  e(t)  and  k|j(t ) ,  Fm=0.1,  Fahl=0.0l,  Fah2=0  004, 
Fah3=0.002,  ^33=0.05,  433=0  05,  4l3=0.0125  and  4ii=0.01,  ^=0.5, 
bh=bc,  ^,=0,  i=l,  2,  £1/£2=2,  £i/£3=4  and  three  £)=£  values;  a)ygia 
versus  fth,  b)  pa  versus 


4.6.  EXPERIMENTAL  VALIDATION 
4.6.1.  Spur  (iear  Pair  Dynamics 

First  we  compare  our  theory  with  experimental  results  of  Kubo  [28],  as  extracted 
from  Reference  [24].  Kubo  used  a  heavily  damped  (£33=0. 1 )  four-square  spur  gear 
test  rig,  and  measured  dynamic  factor  as  the  ratio  of  the  dynamic  to  static  tooth  root 
stresses.  The  experimental  set-up  was  designed  to  support  a  gear  pair  with  very  stiff 
shafts  and  bearings.  Therefore  our  gear  pair  model  can  represent  the  test  rig 
adequately,  similar  to  Chapter  11.  The  static  transmission  error  e(t)  and  tune-varymg 
mesh  stiffness  of  the  tested  gear  pair  has  been  predicted  using  an  existing  elastic 
spur  gear  model  [2 1  ]  and  then  equation  (4.4)  is  solved  to  predict  the  dynamic  response 
Here,  we  define  the  dynamic  factor  as  the  dynamic  to  static  mesh  force  ratio  which  is 
equivalent  to  the  dynamic  facloi  calculation  based  on  the  stress  analysis  under  the 
assumption  that  the  change  in  the  moment  ami  due  to  changes  in  the  contact  point  is 
negligible.  Figure  4.15  compares  the  envelope  of  measurements  obtained  by  Kubo  for 
several  tooth  pairs  with  our  predictions.  When  the  mesh  stiffness  kh(t)  is  assumed  to 
be  tune-invariant,  the  predicted  jump  discontinuity  is  not  as  large  as  the  jump  seen  in 
experimental  data  and  the  predicted  transition  frequency  is  higher  than  the  measured 
value,  as  reported  in  Chapter  II.  However  predictions  improve  significantly  when 
sinusoidally  varying  kh(t)  and  eft )  at  f2h  are  used,  as  given  by  equation  (4.5).  A 
sharp  jump  discontinuity  is  found  which  matches  well  with  experiment,  and  the 
predicted  dynamic  factor  is  very  close  to  measured  envelope.  Finally,  prediction  agrees 
with  experiment  even  better  when  periodically  varying  e(t)  and  kh(t)  with  the  first  three 
Fourier  coefficients  considered.  This  figure  clearly  shows  that  the  time-invariant  linear 
or  non  linear  model  can  not  predict  the  true  dynamic  behavior,  as  time-varying  mesh 
stiffness  must  be  included  in  the  non-linear  mathematical  formulation. 
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Figure  4.15  Comparison  of  theory  with  Kubo's  124,28]  experimental  results. 
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4.6.2.  Geared  Rotor-Bearing  System 

As  the  second  example  case,  experimental  result*  of  Munro  (27]  are  compared 
with  our  geared  rotor-bearing  model  of  Figure  4.2a  and  equation  (4  3).  These 
experimental  results  were  used  earlier  to  validate  time-invariant  multi-degree  of  freedom 
non-linear  system  of  Chapter  III  in  which  we  had  to  reduce  the  damping  ratio  and 
increase  the  excitation  S(t)  in  order  to  correlate  theory  with  experiment.  Now  we  go 
back  to  the  original  system  parameters  given  in  References  (27J  and  Chapter  111,  and 
predict  the  dynamic  transmission  error  spectra  using  both  time-varying  and  lime- 
invariant  gear  mesh  stiffness  formulations  for  a  three  degree  of  freedom  non-linear 
model  with  gear  backlash  and  linear  bearings.  Figure  4,16  compares  results  at  the 
design  load  which  corresponds  to  the  minimum  excitation  f(i ),  In  this  case,  the  mean 
load  to  alternating  load  ratio  F  **  Fm  /  is  very  large,  say  F  «■  30.  Time'  invariant 

stiffness  model  prediction  differs  considerably  from  the  experiment  in  Figure  4,16, 
Predicted  amplitudes  are  considerably  lower  than  the  measurements,  and  a  significant 
jump  discontinuity  found  experimentally  around  the  second  natural  frequency  0)y  Is  not 
even  predicted  by  this  formulation.  But  when  the  periodically-varying  mesh  stiffness  is 
considered,  our  model  predicts  the  frequency  response  accurately  including  the  Jump 
discontinuity.  The  reason  for  a  large  jump  around  cojj  for  a  very  heavily  loaded  sysiem 
(F»  30)  is  now  clear  The  second  natural  frequency  (Ojj  of  the  corresponding  LTI 
system  is  nearly  twice  the  first  natural  frequency  Wj,  This  fortes  the  second  primary 
resonance  at  flj,  »  u)jj  to  coincide  with  the  parametric  resonance  at  Uj(»2cj)j, 

consequently  a  very  large  jump  discontinuity  is  develo|>ed. 

At  3/4  of  the  design  load  with  F  *  10,  the  effect  of  k|,(t)  is  observed  in  l  iguic 
4  17.  Although  the  time-invariant  model  predicts  the  Jump  at  u>j|,  predicted  amplitudes 


Fl|uft  4.16  Coinjwlion  of  theory  with  Murvo'a  (27]  experiment  at  the  dealgn  load. 
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Figure  4. 17  Comparison  of  theory  with  Munro's  [27]  experiment  at  3/4  of  the 
design  load. 
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Figure  4.18  Comparison  of  theory  with  Munro's  [27]  experiment  at  1/2  of  the 
design  load, 
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Figure  4,19  Compariiion  of  theory  with  Munro’s  (27)  experiment  at  1/4  of  the 
design  load. 


2(H) 


are  far  below  the  measurements.  The  inclusion  of  kb(t)  improves  predictions 
drastically.  Similarly  at  1/2  and  1/4  design  loads,  mesh  stiffness  variation  as  shown  in 
Figures  4.18  and  4.19  respectively,  affects  the  frequency  response  significantly  and 
yields  predictions  closer  to  the  measurements. 

4.7.  CONCLUDING  REMARKS 

This  chapter  on  the  non-linear  dynamics  of  a  geared  rotor-bearing  system  with 
time-varying  mesh  stiffness  kb(t),  as  excited  by  the  static  transmission  error  under  a 
mean  load,  has  resolved  a  number  of  fundamental  issues.  First,  the  interaction  between 
time-varying  mesh  stiffness  kh(t)  and  mean  torque  load  has  been  understood. 
Second,  frequency  response  of  the  corresponding  LTV  system  has  been  studied,  and 
the  resonances  associated  with  parametric  and  forced  excitations  have  been  identified. 
Third,  dynamic  interactions  between  kh(t)  and  system  non-linearities  associated  with 
gear  backlash  and  radial  clearances  in  rolling  element  bearings  have  been  investigated;  a 
strong  interaction  between  kh(t)  and  gear  backlash  is  found  where  as  the  coupling 
between  kb(t)  and  bearing  non-linearities  is  weak.  Finally,  our  time-varying  non¬ 
linear  formulation  has  yielded  good  predictions  as  compared  with  benchmark  experiments. 
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CHAPTER  V 


CONCLUSION 


5.1.  SUMMARY 

In  this  study,  linear  and  non-linear  mathematical  models  of  a  generic  geared 

rotor-bearing  system  shown  in  Figure  1  la  are  developed,  and  several  modeling  issues 
which  have  never  been  addressed  previously  in  the  literature  are  investigated  in  depth. 
In  Chapter  I,  a  dynamic  finite  element  model  of  the  system  is  developed.  Effects  of 
several  system  parameters  such  as  torsional  and  transverse  flexibilities  of  the  shafts  and 
prime  mover/load  inertias  are  investigated,  and  modes  of  interest  arc  identified.  Three 
reduced  order  linear  time-invariant  models  are  developed  and  the  conditions  under 
which  such  models  are  suitable  are  determined  by  comparing  the  eigen-solutions  with 
the  finite  clement  model  results. 

In  Chapter  II,  non-linear  frequency  response  characteristics  of  a  spur  gear  pair 
with  backlash  and  time-invariant  mesh  stiffness  are  examined  for  both  external  and 
internal  excitations.  The  internal  excitation  is  of  importance  from  the  high  frequency 
noise  and  vibration  control  view  point  and  it  represents  the  overall  kinematic  or  static 
transmission  error  e(t).  Such  problems  may  be  significantly  different  from  the  rattle 
problems  associated  with  external,  low  frequency  torque  excitation.  Two  solution 
methods,  namely  the  digital  simulation  technique  and  the  method  of  harmonic  balance 
have  been  used  to  develop  the  steady  state  solutions  for  the  internal  sinusoidal 
excitation.  Difficulties  associated  with  the  determination  of  the  multiple  solutions  at  a 
given  frequency  in  the  digital  simulation  technique  have  been  resolved  as  one  must 


search  the  entire  initial  conditions  map.  Such  solutions  and  the  transition  frequencies 
for  various  impact  situations  are  found  analytically  by  the  method  of  harmonic  balance. 
Further,  the  principle  of  superposition  can  be  employed  to  analyze  the  periodic 
transmission  error  excitation  and/or  combined  excitation  problems  provided  the 
excitation  frequencies  are  sufficiently  far  apart  from  each  other.  Predictions  compare 
reasonably  well  with  the  experimental  data  available  in  the  literature. 

In  Chapter  IH,  non-linear  frequency  response  characteristics  of  a  geared  rotor¬ 
bearing  system  are  exainined.  A  three  degree  of  freedom  dynamic  model  is  developed 
which  includes  non-linearities  associated  with  radial  clearances  in  the  radial  rolling 
element  bearings  and  backlash  between  a  spur  gear  pair;  linear  time-invariant  gear 
meshing  stiffness  is  assumed.  Bearing  non-linear  stiffness  function  is  approximated 
for  convenience  sake  by  a  simple  model  which  is  identical  to  that  used  for  the  gear 
mesh.  This  approximate  bearing  model  has  been  verified  by  comparing  the  steady  state 
frequency  spectra.  Applicability  of  both  analytical  and  numerical  solution  techniques  to 
the  multi-degree  of  freedom  non-linear  problem  is  investigated.  Proposed  theory  is 
validated  by  comparing  the  results  with  available  experimental  data.  S'-veral  key  issues 
such  as  non-linear  modal  interactions  and  differences  between  internal  static 
transmission  error  excitation  and  external  torque  excitation  are  discussed.  Additionally, 
parametric  studies  are  performed  to  understand  the  effect  of  system  parameters  such  as 
bearing  stiffness  to  gear  mesh  stiffness  ratio,  alternating  to  mean  force  ratio  and  radial 
bearing  preload  to  mean  force  ratio  on  the  non-linear  dynamic  behavior.  A  criterion 
used  to  classify  the  steady  state  solutions  is  presented,  and  the  conditions  for  chaotic, 
quasi-periodic  and  subharmonic  steady  state  solutions  are  determined.  Two  typical 
routes  to  chaos  observed  in  this  geared  system  are  also  identified. 
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Non-linear  frequency  response  characteristics  of  a  geared  rotor-bearing  system 
with  time-varying  mesh  stiffness  are  examined  in  Chapter  IV.  The  spur  gear  pair 
model  of  Chapters  II  and  geared  rotor-bearing  system  model  of  Chapter  ID  are 
modified  to  include  periodic  mesh  stiffness  k^ft).  Governing  non-linear  time-varying 
equations  which  include  clearance  non-linearities  associated  with  gear  backlash  and 
rolling  element  bearings,  as  excited  by  the  static  transmission  error  e(t)  under  a  mean 
torque  load,  are  solved  using  digital  simulation  technique.  Resonances  of  the 
corresponding  linear  time-varying  (LTV)  system  associated  with  parametric  and 
external  excitations  are  identified  using  the  method  of  multiple  scales  and  digital 
simulation.  Interactions  between  mesh  stiffness  variation  and  clearance  non-linearities 
have  been  investigated;  a  strong  interaction  between  time  varying  mesh  stiffness  kjj(t) 
and  gear  backlash  is  found  whereas  the  coupling  between  kh(t)  and  bearing  non- 
linearities  is  weak.  The  predictions  yielded  by  the  proposed  time-varying  non-linear 
model  agree  well  with  the  experimental  results  available  in  the  literature. 

5.2  FUTURE  RESEARCH  AREAS 

The  following  topics  are  identified  as  areas  of  future  research,  based  on  the 
present  study  on  non-linear  dynamic  analysis  of  geared  rotor-bearing  systems: 

1 .  Extension  of  current  models  to  study  dynamics  of  helical  and  bevel  gear  drives, 
and  multi-gear  mesh  systems  such  as  planetary  gear  trains. 

2.  Iteration  procedures  for  the  transmission  error  input  and  inclusion  of  the 
interactions  between  external  and  internal  kinematic  error  excitations.  Also 
inclusion  of  side-bands  in  non-linear  dynamic  analysis. 
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3 .  Design  of  a  passive  vibration  control  strategy  using  squeeze  film  dampers.  And, 
investigation  of  semi-active  and  active  vibration  and  noise  control  schemes  using 
piezoelectric  actuators  or  hydraulic  mounts. 

4 .  Statistical  energy  analysis  of  non-linear  geared  systems. 
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APPENDIX  A 


USER'S  GUIDE  FOR  THE  GEARED  ROTOR  DYNAMICS 

PROGRAM  -  GRD 

A.l.  DESCRIPTION 

Geared  Rotor  Dynamics  Program  "GRD"  is  a  general  purpose  finite  element 
computer  program  to  analyze  dynamics  of  a  system  consisting  of  two  shafts  supported 
on  bearings  and  coupled  by  a  gear  mesh.  It  computes  natural  frequencies,  mode  shapes 
and  the  vibration  response  of  the  system  to  static  transmission  error  excitation,  mass 
unbalances  and  geometric  eccentricities  of  gears.  Many  rigid  disks,  bearings,  and 
hollow  shafts  can  be  considered  by  GRD.  Since  GRD  uses  finite  elements  method,  the 
system  should  first  be  discretized  to  small  rotor  elements.  It  is  also  necessary  to  enter 
data  for  each  element  in  sequence. 

A. 2.  USING  THE  PROGRAM 

There  are  two  ways  to  input  data  in  GRD.  Data  can  be  entered  through  a  file  or 
interactively.  The  first  set  of  questions  will  give  the  user  this  option  (the  Primary 
Option).  Display  on  screen. 

ENTER  DATA  INTERACTIVELY  OR  FROM  A  DATA  FILE? 

OPTION  (1) :  ENTER  DATA  INTERACTIVELY 

OPTION  (2) :  ENTER  DATA  FROM  FILE 

CHOOSE  OPTION  1  OR  2  (NO  DEFAULT)  »  _ 

A. 2.1.  OPTION  I:  Interactive  Data  Entry 

By  selecting  this  option,  the  user  will  be  prompted  to  enter  all  the  data  that  will  be 
necessary  to  nm  this  program.  In  the  process,  a  new  data  file  will  be  created  from  the 
user's  input.  The  data  entered  will  be  reorganized  in  a  format  so  that  the  program  may 
read  from  it.  For  this  reason,  a  new  data  file  will  be  created  and  the  user  will  be 
prompted  to  name  this  new  file  as  follows: 

DATA  WILL  BE  INPUT  INTERACTIVELY,  AND  IN  THE  PROCESS.  YOUR  INPUT 
WILL  BE  ORGANIZED  IN  A  NEW  DATA  FILE  FOR  THIS  PROGRAM  TO 
READ  FROM. 

ENTER  AN  INPUT  FILE  NAME  »  _ 
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If  an  error  should  arise  from  opening  this  new  file,  the  user  will  be  given  the  option  to 
try  again: 

♦ERROR  IN  OPENING  FILE*  DO  YOU  WISH  TO  TRY  AGAIN? 

TYPE  (1)  FOR  YES  OR  (2)  FOR  NO  »_ 

If  YES,  then  the  user  will  be  asked  the  file  name  again.  If  NO,  the  program  will 
return  to  the  Primary  Option. 


2.1.1.  Data  Entry  (Groupings 

In  this  section,  the  user  will  enter  the  necessary  data.  The  questions  will  be 
grouped  into  the  following  categories. 


1.  GENERAL  DATA 

2.  MATERIAL  PROPERTIES 

3.  ELEMENT  PROPERTIES 

4.  GEAR  MESH  PROPERTIES 

5.  FORCED  RESPONSE  DATA  (If  user  chooses  this  option) 

J.  General  Data: 

In  this  section,  the  following  general  information  about  the  overall  system  will  be 
entered: 


1)  THE  NUMBER  OF  ROTOR  ELEMENTS  IN  THE  FIRST  SHAFT 

2)  THE  NUMBER  OF  ROTOR  ELEMENTS  IN  THE  SECOND  SHAFT 

3)  THE  NUMBER  OF  ROTOR  ELEMENTS  BEFORE  THE  OEAR  IN  THE  FIRST  SHAFT 

4)  THE  NUMBER  OF  ROTOR  ELEMENTS  BEFORE  THE  OEAR  IN  THE  SECOND 
SHAFT 


As  an  example,  consider  the  system  shown  in  Figure  Al.  The  fust  step  is  to  the 
divide  rotors  into  small  pieces  (finite  rotor  elements)  to  obtain  a  finite  element  model  of 
the  system  (as  illustrated  in  Figure  A2).  The  shorter  the  rotor  elements,  the  higher  the 
accuracy  and  the  longer  the  computation  time  In  Figure  A2,  each  shaft  is  divided  into 
4  pieces.  The  numbering  can  start  at  any  end  of  the  driving  rotor  and  finish  at  any  end 
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PINION 


Figure  Al.  A  Typical  Geared  Rotor  System 
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Figure  A2.  Finite  Element  Model 
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of  the  driven  rotor.  In  Figure  A2  for  example,  we  start  at  the  left  end  of  the  first  shaft 
and  finished  numbering  at  the  right  end  of  the  second  shaft.  According  to  this 
configuration  (Figure  A2),  the  number  of  rotor  elements  in  the  first  and  second  shafts 
are  both  4  and  the  number  of  rotor  elements  before  the  pinion  and  gear  are  2  and  6, 
respectively. 


2.  Material  Progenies: 

In  this  section,  material  properties  of  the  shaft  considered  will  be  entered.  The 
following  will  be  required: 


1 )  THE  DENSITY  OF  THE  SHAFT  MATERIAL  (kg/m*) 

2)  THE  VISCOUS  DAMPING  COEFFICIENT  (s) 

3)  THE  HYSTERET1C  DAMPING  COEFFICIENT 

4)  THE  MODULUS  OF  ELASTICITY  (N/m2) 

5)  SHEAR  MODULUS  OF  ELASTICITY  (N/m2) 


3.  Element  Eigpcciga: 

The  main  portion  of  data  to  be  entered  is  in  tikis  section, 
following: 

1)  DEFINING  THE  DISK  ELEMENT 

2)  DEFINING  THE  BEARING  ELEMENT 

3)  DEFINING  THE  ROTOR  ELEMENT  (also  the  default  selection) 

4)  DEFINING  THE  GEAR  ELEMENT 

5)  FINISH  WITH  DATA  INPUT  FOR  FIRST  SHAFT 

6)  FINISH  WITH  DATA  INPUT  FOR  SECOND  SHAFT 


(DEFAULT=7800) 

(DEF=0) 

(DEF=0) 

(DEF=.20?E12) 

(DEF=.0795E12) 


The  options  given  are  the 


The  most  important  part  in  this  section  of  questions  is  the  order  of  data  entry.  The 
order  is  based  on  the  numbered  finite  element  model  (see  Figure  A2).  For  the  finite 
element  model  shown  in  Figure  A2,  the  order  of  data  input  is  as  follows: 

1)  DEFINING  BEARING  #1 

2)  DEFINING  ROTOR  #1 

3)  DEFINING  ROTOR  #2 

4)  DEFINING  GEAR#1 

5)  DEFINING  ROTOR  #3 

6)  DEFINING  ROTOR  #4 

7)  DEFINING  BEARING  #2 

8)  TYPE  IN  OPTION  5  (FINISHED  WITH  DATA  INPUT  FOR  SHAFT  #1 ) 

9)  DEFINING  BEARING  #3 

10)  DEFINING  ROTOR  #5 
1DDEHNING  ROTOR  #6 
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12)  DEFINING  GEAR#2 

13)  DEFINING  ROTOR  #7 

14)  DEFINING  ROTOR  #8 

15)  DEFINING  BEARING  #4 

16)  TYPE  IN  OPTION  6  (FINISHED  WITH  DATA  INPUT  FOR  SHAFT  #2) 


From  the  example  above,  one  can  see  ihe  ordering  of  data  entry  moving  from  the  First 
element  in  shaft  #1  to  the  last  element  in  shaft  #2.  Ultimately,  the  program  will  also 
read  from  the  data  file  in  this  order. 

According  to  the  option  selected  (disk,  gear,  rotor  element  or  bearing),  the 
following  information  is  required: 


The  disk  element  -  option  1: 


1 )  DISK  OUTER  DIAMETER  IN  (m) 

2)  WIDTH  OR  THICKNiiSS  OF  DISK  IN  (m) 

3)  MATERIAL  DENSITY  OF  DISK  IN  (kg/m3)  (DEF=7800) 


The  bearing  element  -  option  2: 


THE  STIFFNESS  (N/mt 

1)  K*x 

2)  K*y 

3) Kyx 

4)  Kyy 

THE  DAMPING  (N-s/m) 

5) 

6)  CXy 

7)  CyX 

8)  Cyy 


The  rotor  element  -  option  3: 


1 )  LENGHT  OF  ROTOR  ELEMENT  IN  (m) 

2)  OUTER  DIAMETER  OF  ROTOR  ELEMENT  IN  (m) 

3)  INNER  DIAMETER  OF  ROTOR  ELEMENT  IN  (m)  (for  hollow  rotor)  (DEF«0) 

4)  THE  AXIAL  LOAD  IN  (N)  (DEF=0) 


The  gear  element  -  option  4: 

J )  PITCH  CIRCLE  DIAMETER  IN  (m) 
2)  FACE  WIDTH  IN  (m) 
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3)  GEAR  MATERIAL  DENSITY  IN  (kg/m3)  (DEF-7800) 


Option  5  should  be  chosen  when  the  elements  of  the  first  shaft  have  been  defined. 
The  user  will  then  be  given  the  material  properties  of  the  second  shaft.,  and  given  the 
option  to  make  changes  to  this  set  of  data.  Similarly,  option  6  should  be  entered 
when  all  elements  of  the  second  shaft  have  been  entered. 


4.  Gear  Mesh  Properties: 

The  following  gear  mesh  properties  will  be  entered  in  this  section: 


1)  GEAR  MESH  STIFFNESS  (N/m) 

2)  GEAR  MESH  DAMPING  (N-s/m)  (DEF=0) 

3)  BASE  CIRCLE  DIAMETER  OF  FIRST  GEAR  (m) 

4)  BASE  CIRCLE  DIAMETER  OF  SECOND  GEAR  (m) 


Following  the  input  of  gear  mesh  data,  the  program  will  also  ask  the  user  for  the 
number  of  natural  frequencies  wanted  in  the  output. 


5.  Forced  Response  Data; 

Here,  the  user  is  given  the  option  to  have  forced  response  calcula»ed.  If  the 
selection  is  NO,  then  interactive  input  will  be  complete.  If  the  selection  is  YES,  then 
the  following  options  are  given: 


1 )  WHIRLING  ORBIT  AT  A  SPECIFIED  NODE 

2)  DEFLECTIONS  AT  A  SPECIFIED  NODE 

3)  DYNAMIC  LOAD  TO  STATIC  LOAD  RATIO  AT  THE  MESHING  POINT 


Following  the  options  menu,  the  general  information  for  this  section  will  be  required: 


1 )  GEOMETRIC  ECCENTRICITY  (RUNOUT)  OF  GEAR  I  (m)  (DEF-O) 

2)  GEOMETRIC  ECCENTRICITY  (RUNOUT)  OF  GEAR  2  (m)  (DEF-O) 

3)  PEAK  TO  PEAK  VALUE  OF  STATIC  TRANSMISSION  ERROR  (m) 

4)  MASS  UNBALANCE  OF  GEAR  1  (kg-m)  (DEF-O) 

5)  MASS  UNBALANCE  OF  GEAR  2  (kg-m)  (DEF-O) 

6)  NUMBER  OF  TEETH  IN  GEAR  1  (PINION) 

7)  AVERAGE  FORCE  TRANSMITTED  (N) 

8)  VALUE  OF  MODAL  DAMPING  (DEF=.0 1 ) 


After  the  general  information  for  this  section  has  been  completed,  the  program  will 
require  certain  information  based  from  the  choice  selected  from  the  options  menu  : 
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Forced  response  option  1:  Whirling  orbit  at  a  specified  node,  the  following 
information  will  be  required. 

1 )  NUMBER  OF  NODE  AT  WHICH  WHIRL  ORBIT  IS  REQUIRED 

2)  ROTATIONAL  SPEED  OF  SHAFT  1  (rad/s) 


Forced  response  option  2:  Deflections  at  a  specified  node  will  be  required. 


1)  STARTING  ROTATIONAL  SPEED  OF  SHAFT  I  (rad/s)  (DEF=0) 

2)  UPPER  LIMIT  OF  ROTATIONAL  SPEED  OF  SHAFT  1  (rad/s) 

3)  INCREMENT  FOR  ROTATIONAL  SPEED  (rad/s)  (DEF=5) 

4)  DIRECTION  AT  WHICH  FORCED  RESPONSE  IS  WANTED  ( 1 ,2  OR  3)  (DEF=1) 

5)  NUMBER  OF  NODE  AT  WHICH  FORCED  RESPONSE  IS  WANTED 


In  question  (4),  enter:  I  for  deflection  in  pressure  line  direction;  2  for  deflection  in 
direction  perpendicular  to  pressure  line;  3  for  torsional  deflections 

1)  STARTING  ROTATIONAL  SPEED  OF  SHAFT  1  (rad/s)  (DEF=0) 

2)  UPPER  LIMIT  OF  ROTATIONAL  SPEED  OF  SHAFT  1  (rad/s) 

3)  INCREMENT  FOR  ROTATIONAL  SPEED  (rad/s)  (DEF=5) 

Having  completed  this  section  on  forced  response,  the  program  will  then  read  from  the 
new  data  file  created  from  the  interactive  input.  The  output  will  be  in  an  output  data  file 
namedFOR001.DAT  . 

A.2.2.  OPTION  II:  Enter  Data  From  File 

With  this  selection,  the  user  is  then  requested  to  give  the  name  of  the  file  where  the 
data  is  stored,  instead  of  entering  data  interactively.  The  program  will  then  proceed  to 
read  from  the  file.  If  an  error  should  occur  in  opening  the  file,  then  the  user  will  be 
given  the  option  to  try  again.  If  the  user  opts  to  try  again,  then  the  user  will  be  asked 
the  file  name  again.  If  no,  then  the  program  will  return  to  the  primary  option.  After  the 
program  has  read  from  the  file  the  output  will  be  found  in  an  output  data  file. 
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A. 4.  SAMPLE  OUTPUT  FILE  FOROOJ.DAT 


GEARED  ROTOR  DYNAMICS  PROGRAM 
GRD 

GEAR  DYNAMICS  AND  GEAR  NOISE  RESEARCH  LABORATORY 
THE  OHIO  STATE  UNIVERSITY 

ft******************************************************************** 


THE  FIRST  SHAFT: 


MATERIAL  PROPERTIES  OF  THE  SHAFT: 


DENSITY  OF  THE  MATERIAL . « 

ELASTIC  MODULUS . - 

SHEAR  MODULUS . - 

VISCOUS  DAMPING  COEFFICIENT . - 
HYSTERETIC  LOSS  FACTOR . - 


0 . 78000E+04  KG/M*  *  3 
0.20700K+12  N/M* * 2 
0 . 79500E+11  N/M*  *2 
0 . 00000E+00  S 
0 . 00000E+00 


AT  Z-  0.0000  M  THERE  EXISTS  A  BEARING  WITH  THE 
FOLLOWING  STIFFNESS  AND  DAMPING  COEFFICIENTS: 


KXX- 

0 . 10000E+10 

N/M 

KXY- 

o 

o 

+ 

u 

o 

o 

o 

o 

o 

o 

N/M 

KYX- 

O.OOOOOE+OO 

N/M 

KYY- 

0. 10000E+10 

N/M 

CXX- 

O.OOOOOE+OO 

N-S/M 

CXY- 

O.OOOOOE+OO 

N-S/M 

CYX- 

0  .  OOOOCE+OO 

N-S/M 

CYY«= 

o 

O 

+ 

U 

o 

o 

o 

o 

o 

o 

N-S/M 

AT  Z-  0.0000  M  THERE  EXISTS  A  FINITE  ROTOR  ELEMENT  WITH  THE 
FOLLOWING  SPECIFICATIONS: 


LENGTH  OF  THE  ELEMENT 

OUTER  DIAMETER . 

INNER  DIAMETER . 

AXIAL  LOAD . 


0 . 50000F.-01  M 
0.30000E-01  M 
0 . 00Q00E+00  M 
0 . 0000OE+00  N 


AT  Z-  0.0500  M  THERE  EXISTS  A  riNITF.  ROTOR  ELEMENT  WITH  THE 
FOLLOWING  SPECIFICATIONS: 


LENGTH  OF  THE  ELEMENT 

OUTER  DIAMETER . 

INNER  DIAMETER . 

AXIAL  LOAD . 


0.50000E-01  M 
0 . 30000E-0I  M 
0 . 00000E+00  M 
O.OOOOOE+OO  N 


AT  Z-  0.1000  M  THERE  EXISTS  A  RIGID  DISK  WITH  THE  FOLLOWING 
SPECIFICATIONS: 


OUTER  DIAMETER . -  0.13500E+00  M 
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WIDTH . -  0.25400E- 01  H 

MATERIAL  DENSITY.  ,  ,  0.78000E+04  KO/M»*3 


AT  E-  0.1000  M  THERE  EXISTS  A  flHITE  ROTOR  EL&R05HT  WITH  THE 
POLLOMXNG  SPECIFICATIONS t 


LENGTH  Or  THE  ELEMENT 

OUTER  DIAMETER . 

INNER  DIAMETER . 

AXIAL  LOAD . 


0 . 50000E-01  M 
0.30000E-0)  M 
0. 00000E+00  M 
0 . 00000E+00  N 


AT  Z-  0,1500  M  THERE  EXISTS  A  FINITE  ROTOR  ELEMENT  WITH  THE 
PQLLOWINU  SPECiriCATIONS: 


LENOTH  Or  THE  ELEMENT 

OUTER  DIAMETER . 

INNER  DIAMETER . 

AXIAL  LOAD . 


0 , 50000E-01  M 
0 , JOO&OE-Ol  M 
0 , OOOOOE+  00  M 
0.00000E400  N 


AT  E-  0,2000  M  THERE  EXISTS  A  EEAKINO  WITH  THE 
FOLLOWING  STIFFNESS  AND  DAMPING  COEFFICIENTS i 


KXX- 

0. 1 OOOOR4 10  N/M 

KXY- 

0 , OOOOOK  4  00 

N/M 

KYX- 

0 ,  OOOOOE  +  OO  N/M 

KYY- 

0. JO00OE+10 

N/M 

CXX- 

0. OOOOOE  +  OO  N-S/M 

CXY- 

0. O0000E+00 

N-S/M 

CYX- 

C  .  00000E4  00  N-S/M 

CYY- 

O  ,  OOOOO E- +  00 

II-*/ M 

THE  SECOND  SHAFT: 


MATERIAL  PROPERTIES  OP  THE  SHAFT: 


DENSITY  Of  THE  MATERIAL . - 

ELASTIC  MODULUS . - 

SHEAR  MODULUS.  .  ,  . . - 

VISCOUS  DAMPING  COEFFICIENT  ,  - 
HY8TERETIC  LOSS  PACTOR . - 


0,7B000F,  +  04  KG/M*  *  3 
0.2070QE412  N/M* * 2 
0.79S00E+1)  N/M* * 2 
0 , 00000E4  00  8 
0 , 00QQ0E+  00 


AT  Z-  0,0000  M  THERE  EXISTS  A  REARING  WITH  THE 
POLLOWINO  STIFFNESS  AND  DAMPING  COEfPlCIENTBl 


KXX» 

0.10000E+10 

N/M 

KXY* 

0. OOOOOE+OO 

N/M 

KYX- 

0. 00000E  +  00 

N/M 

KYY- 

0, 10000E+10 

N/M 

CXX- 

0 .  OOOOOE  +  OO 

N-S/M 

CXY- 

0, OOOOOE+OO 

N-S/M 

CYX- 

0, OOOOOE+OO 

N-S/M 

CYY- 

0, OOOOOE+OO 

N-S/M 

AT  Z-  0,0000  M  THERE  EXISTS  A  PINITE  ROTOR  ELEMENT  WITH  THE 
POLLOW1NG  6PEC7 PICATIONS I 


LENGTH  OP  THE  ELEMENT 

OUTER  DIAME'ltR . 

INNER  DIAMETER . 

AXIAL  LOAD . 


0.50000E-01  M 
0 . 40000P.-01  M 
0 i OOOOOE+OO  M 
0, OOOOOE+OO  N 
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AT  Z-  O.OSOO  M  THERE  EXISTS  A  FINITE  ROTOR  ELEMENT  WITH  THE 
FOLLOWING  SPECIFICATIONS : 


LENGTH  Or  THF  ELEMENT 

OUTER  DIAMETER . 

INNER  DIAMETER . 

AXIAL  LOAD . 


0.50000E-01  M 
0.40000E-01  M 
O.OOOOOE+OO  M 
0 .  OOOOOE+OO  N 


AT  Z-  0.I0C0  M  THERE  EXISTS  A  RIGID  DISK  WITH  THE  FOLLOWING 
SPECIFICATIONS : 


OUTER  DIAMETER . 

WIDTH . 

MATERIAL  DENSITY. . . . 


0.24500E+00  M 
0 . 25400E-01  M 
0 . 78OO0E+O4  KG/M* *3 


AT  Z-  0.1000  M  THERE  EXISTS  A  FINITE  ROTOR  ELEMENT  WITH  THE 
FOLLOWING  SPECIFICATIONS: 


LENOTH  OF  THE  ELEMENT 

OUTER  DIAMETER . 

INNER  DIAMETER . 

AXIAL  LOAD . 


0.50000E-01  M 
0 . 40000E-01  M 
O.OOOOOE+OO  M 
O.OOOOOE+OO  N 


AT  Z-  0.1500  M  THERE  EXISTS  A  TINITE  ROTOR  ELEMENT  WITH  THE 
FOLLOWING  SPECIFICATIONS: 


LENGTH  Or  THE  ELEMENT 

OUTER  DIAMETER . 

INNER  DIAMETER . 

AXIAL  LOAD . 


0 . 50000E-01  M 
0 . 40000E-01  M 
O.OOOOOE+OO  M 
O.OOOOOE+OO  N 


AT  Z-  0.2000  M  THERE  EXISTS  A  BEARING  WITH  THE 
FOLLOWING  STITFNESS  AND  DAMPING  COEFFICIENTS: 


KXX-  0 . 10000E+10 
KYX-  O.OOOOOE+OO 
CXX-  O.OOOOOE+OO 
CYX-  O.OOOOOE+OO 


GEAR  MESH  PROPERTIES: 


N/M 

KXY- 

N/M 

KYY- 

N-S/M 

CXY- 

N-S/M 

CYY- 

O.OOOOOE+OO  N/. 

0 . 10000E+10  N/M 
O.OOOOOE+OO  N-S/M 
O.OOOOOE+OO  N-S/M 


BASE  CIRCLE  DIAMETER  OF  GEAR  I...- 
BASE  CIRCLE  DIAMETER  Or  GEAR  2...- 

AVERAGE  MESH  STIFFNESS . - 

AVERAGE  MESH  DAMPING . - 


0.1270  M 
0.2310  M 
0 . 200C0E+09  N/M 
O.OOOOOE+OO  N-S/M 


a***********************"******************************************** 


A.  TREE  VIBRATION  ANALYSIS: 


FOLLOWING  NATURAL  FREQUENCIES  AND  MODESHAPES  ARE  CALCULATED 
TOR  THE  ABOVE  SPECIFIED  SYSTEM: 
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MODE  1 


NATURAL  FREQUENCY-  0.0002  (HZ.) 

CORRESPONDING  MODESHAPES  : 


in  y 

DISP.  IN  X 

ROT.  ABOUT  Y 

ROT.  ABOUT  X 

ROT.  ABOUT  Z 

- . 260E-08 

0 . 000E+00 

0.  OOOE+OO 

-. 105E-06 

0 . 599E+01 

-.358E-08 

0 . 000E+00 

0. OOOE+OO 

0.108E-06 

0 . 599E+01 

0 . 119E-08 

0 . OOOE+OO 

0. OOOE+OO 

-.411E-07 

0 . 599E+01 

- . 524E-08 

0 .  OOOE+OO 

0. OOOE+OO 

-. 983E-07 

0 . 599E+01 

- . 299E-08 

0 . 000E+00 

0. OOOE+OO 

0 . 136E-06 

0 . 599E+01 

0 . 771E-08 

0  .OOOE+OO 

0 .OOOE+OO 

0 . 499E-06 

0.329E+01 

0 . 253E-07 

0 . OOOE+OO 

0. OOOE+OO 

0 . 53SE-07 

0 . 32  9E+01 

-. 188E-09 

0. 000E+00 

0 .OOOE+OO 

-• .  1 1 8E-05 

0 . 329E+01 

-.256E-07 

0  .  OOOE+OO 

0. OOOE+OO 

0.557E-07 

0 . 329E+01 

- . 781E-08 

0 . OOOE+OO 

0. OOOE+OO 

0 . 505E-06 

0. 329E+01 

A******************************************************************** 

MODE  2 


NATURAL  FREQUENCY-  486.0262  (HZ.) 
CORRESPONDING  MODESHAPES: 


DISP.  IN  Y 

DISP.  IN  X 

ROT.  ABOUT  Y 

ROT.  ABOUT  X 

ROT.  ABOUT  Z 

0. 929E-02 

0 . OOOE+OO 

0. OOOE+OO 

0. 551E+01 

- . 646E+01 

0 . 261E+00 

0  .OOOE+OO 

0. OOOE+OO 

0. 410E+01 

-. 645E+01 

0 . 375E+00 

0 .OOOE+OO 

0. OOOE+OO 

-. 119E-06 

643E+01 

0 . 261E+00 

0  .  OOOE+OO 

0. OOOE+OO 

- . 410E+01 

-. 645E+01 

0 . 929E-02 

0  .  OOOE+OO 

0. OOOE+OO 

- . 551E+01 

- . 646E+01 

-. 950E-02 

0  .  OOOE+OO 

0. OOOE+OO 

-. 180E+01 

0. 110E+01 

-. 918E-01 

0  .OOOE+OO 

0. OOOE+OO 

- . 134E+01 

0 . 110E+01 

- . 129E+00 

0 . OOOE+OO 

0. OOOE+OO 

0 . 995E-09 

0 . 109E+01 

-. 918E-01 

0. OOOE+OO 

0. OOOE+OO 

0 . 134E+01 

0. 110E+01 

- . 950E-02 

0. OOOE+OO 

0. OOOE+OO 

0.180E+01 

0. 110E+01 

MODE  10 


NATURAL  FREQUENCY™  2479.9504  (HZ.) 
CORRESPONDING  MODESHAPES: 


DISP.  IN  Y 

DISP.  IN  X 

ROT.  ABOUT  Y 

ROT.  ABOUT  X 

ROT.  ABOUT  Z 

0 . 124E-01 

0 . OOOE+OO 

0. OOOE+OO 

0 . 459E+01 

0 . 978E+01 

0 . 212E+00 

O.OOOEt-OO 

0 . OOOE+OO 

0 . 288E+01 

0 . 950E+01 

0 . 279E+00 

0. OOOE+OO 

0 . OOOE+OO 

- . 897E-10 

0 . 865E+01 

0 . 212E+00 

0. OOOE+OO 

0. OOOE+OO 

- . 288E+01 

0 . 950E+01 

0. 124E-01 

0. OOOE+OO 

0. OOOE+OO 

~. 459E+01 

0 . 978E+01 

-. 990E-02 

0. OOOE+OO 

0. OOOE+OO 

-.134F.+01 

- . 1 66E+01 

-. 691E-01 

0. OOOE+OO 

0 . OOOE+OO 

- . 897E+00 

- . 162E+01 

- . 917E-01 

0. OOOE+OO 

0. OOOE+OO 

- . 562E-10 

-. 147E+01 

- . 691E-01 

0. OOOE+OO 

0 . OOOE+OO 

0.897E+00 

-. 162E+01 

990E-02 

0. OOOE+OO 

0 . OOOE+OO 

0 . 134E+01 

- . 166E+01 

*******  ************************************************************** 


B.  TORCED  VIBRATION  ANALYSIS: 


FORCED  RESPONSE  AT  THE  SPECIFIED  FREQUENCY  RANGE: 

0 . 591142E-04 
0. 675:86E-04 
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APPENDIX  B 


DIMENSIONAL  AND  NON-DIMENSIONAL  QUANTITIES 
FOR  THE  EXPERIMENTAL  TEST  RIOS 


Table  Bl.  Dimensionless  and  dimensional  parameters  for  Kubo’s  test  rig  [28], 


DIMENSIONAL  QUANTITIES 


Number  of  teeth 

25/25 

Igi.Ig2  [Lg-m2  (lb-ft2)] 

0.00115  (0.0278) 

mci  fkgl 

0.23 

dgi,dg2  [m  (in)] 

0.094  (3.7) 

kb 

rigid 

kh  [N/m  (lbf/in)] 

3.8E8  (2.17E6) 

Fm  [N  (lbDl 

2295  (1030) 

e  [m  (in)] 

1.92E-6  (7.56E-5) 

bfa  [m  (in)] 

0.1  E-3  (0.0039) 

DIMENSIONLESS  QUANTITIES 


c 

(page  44) 

0.1 

Fm 

(page  44) 

0.06 

Fah 

(page  44) 

0.0192 

Fm/Fah 

3.12 
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Table  B2.  Dimensionless  and  dimensional  parameters  for  Munro's  test  rig  [27]. 


DIMENSIONAL  QUANTITIES 


number  of  teeth 
Igl.Ig2  frg-m2  (lb-ft2)] 
mgl,mg2  [kg  (lb)] 
dgi.dgi  [m  (in)] 
bh  [m  (in)] 

32/32 

0.136  (3.29) 

31.1  (68.5) 

0.2  (8) 

0.12E-3  (4.7E-3) 

Design  Load  (DL) 

3/4  DL 

1/2  DL 

1/4  DL 

Fm  [N  (lbf/in  FW)] 

3782  (1700) 

2836(1275)  1891  (850) 

947 (425) 

kh  [N/m] 

3.44E8 

3.22E8 

3.01E8 

2.72E8 

e  [m] 

3.5E-7 

1.06E-6 

1.78E-6 

2.36E-6 

DIMENSIONLESS  QUANTITIES 


Design  Load  (DL) 

3/4  DL 

1/2  DL 

1/4  DL 

Fm  (equation  3.3i) 

0.183 

0.146 

0.105 

0.058 

Pah  (equation  3.3j) 

0.0058 

0.0178 

0.0296 

0.0393 

k11*k22  (equation  3.3 f) 

0.950 

0.966 

0.983 

1.007 

k13-  k23  (equation  3.3g) 

0.242 

0.242 

0.242 

0.242 

C11.C22  (equation  3.3c) 

0.01 

0.01 

0.01 

0.01 

C13,  C23  (equation  3.3d) 

0.00375 

0.00375 

0.00375 

0.00375 

£33  (equation  3.3e) 

0.015 

0.015 

0.015 

0.015 
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APPENDIX  C 


LENGTH  OF  TRANSIENT  SOLUTION  AND  CPU  TIME 

C.l.  Length  of  Transient  Solution 

It  is  necessary  to  run  the  digital  simulation  routines  for  many  cycles  in  order  to 
reach  the  steady  state  solution.  This  can  easily  be  detected  by  a  time  history  repeating 
itself  and  a  phase  plane  trajectory  followed  at  each  period  of  the  steady  state  solution. 
Following  are  found  to  be  the  major  factors  which  determine  the  length  of  the  transient 
region  of  the  time  history: 

1.  Damping  Ratio:  Damping  ratio  is  inversely  proportional  to  the  time  needed  to 
reach  steady  state  solution.  Larger  the  damping,  shorter  the  transients.  Figure  B1 
shows  this  relationship  for  the  single  degree  of  freedom  model  of  gear  pair.  A  similar 
trend  is  observed  in  Multi-degree  of  freedom  model  also. 

2.  Excitation  Frequency:  The  transient  solution  is  found  to  be  longer  in  the 
vicinity  of  the  resonance  frequency  than  the  off-resonance  region.  Figure  B1  compares 
length  of  transient  solution  versus  damping  ratio  curves  corresponding  to  resonance 
and  off-resonance  excitation  frequencies. 

3.  Initial  Conditions:  Transients  are  longer  when  the  initial  displacement  and 
velocity  are  away  from  the  mean  values  of  the  steady  state  solution  sought. 

C.2.  CPU  Times 

CPU  time  needed  in  digital  simulation  depends  on:  i.  the  number  of  increments 
per  period  of  forcing  function  (20-40  points  per  period  is  good  enough),  ii.  the 
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80 


damping  ratio 


Figure  Cl.  Length  of  the  transient  solution  in  a  gear  pair  as  n  function  of  damping 
ratio  for  off-resonance  (fl=0  5)  and  resonance  (Q=0.7)  frequencies 


tolerance  of  the  solution  accuracy  (lxl0‘5  to  IxlO'9  is  acceptable),  and  iii.  number  of 
degrees  of  freedom.  For  the  single  degree  of  freedom  model  with  a  tolerance  of  lxl  O'9 
and  an  increment  of  20  points/period,  nearly  10  sec.  of  CPU  time  is  required  to  run  for 
100  cycles.  The  three  degree  of  freedom  model  with  same  parameters  uses  1.2  min.  of 
CPU  time  for  100  cycles, 
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APPENDIX  I) 


COMPARISON  OF  PREDICTIONS  WITH  NASA  MEASUREMENTS 

In  this  section,  a  comparison  of  the  FEM  predictions  of  the  NASA  Gear  Noise 
Test  Rig  shown  in  Figure  D1  with  the  measurements  of  NASA  is  given.  System 
parameters  and  the  predicted  natural  modes  of  the  system  have  already  been  given  in 
Tables  1.1  and  1.2,  respectively.  The  static  transmission  error  and  mesh  stiffness 
variation  predictions  are  also  shown  in  Figure  4.8.  Here,  the  comparison  of  the 
predicted  frequency  response  spectrum  with  measurements  of  NASA  is  limited  to 
torsional  (angular  acceleration)  vibrations  since  the  set-up  was  not  equipped  for 
transverse  vibration  measurements.  Only  a  qualitative  agreement  is  observed  between 
the  predictions  and  measurements  as  shown  in  Figure  D2  due  to  the  uncertainties 
associated  with  the  sensitivities  of  instruments  and  due  to  the  errors  involved  in 
measurement  such  as  resonance  and  calibration  problems  of  angular  transducers. 


Figure  Dl.  NASA  Gear  Noise  Test  Rig. 
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angular  acceleration  angular  acceleration 


(a) 


0  2000  4000  6000 


frequency  (Hz) 


Figure  D2.  Angular  acceleration  spectrum  of  the  NASA  Gear  Noise  Test  Rig  at  pinion 
location;  a)  measurements  b)  FEM  predictions. 
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APPENDIX  D 


COMPARISON  OF  PREDICTIONS  WITH  NASA  MEASUREMENTS 

In  this  section,  a  comparison  of  the  FEM  predictions  of  the  NASA  Gear  Noise 
Test  Rig  shown  in  Figure  D1  with  the  measurements  by  NASA  is  given.  System 
parameters  and  the  predicted  naturea!  modes  of  the  system  have  already  been  given 
in  Tables  1.1  and  1.2,  respectively.  The  static  transmission  error  and  mesh  stiffness 
variation  predictions  are  also  shown  in  Figure  4.8.  Here,  the  comparison  of  the  predicted 
frequency  response  spectrum  with  measurements  by  NASA  is  limited  to  torsional 
(angular  acceleration)  vibrations.  Only  a  qualitative  agreement  is  observed  between 
the  predictions  and  measurements  as  shown  in  Figure  D2  Quantitative  discrepancies 
are  attributed  to  the  resonance  problems  experienced  with  the  angular  transducers. 


gear  box 


Figure  Dl.  NASA  Gear  Noise  Test  Rig. 
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